Dados geográficos no R
Introdução
Nesta seção, vamos fazer uma breve introdução aos principais conceitos sobre o funcionamento de dados geográficos no R: formatos de dados vetoriais e pacotes; formatos de dados raster e pacotes; Sistemas de Referências de Coordenadas e unidades; fontes de dados geográficos e pacotes; importar e exportar dados geográficos; descrição de objetos espaciais; reprojeção de dados geográficos; e principais operações com dados geográficos. Num segundo momento, iremos criar mapas com seus principais elementos. Por fim, apresentaremos exemplos de aplicações de análises espaciais para dados ecológicos, focadas em resumir informações sobre a biodiverisdade, preparar de dados para compor variáveis preditoras, e como fazer predições espaciais contínuas de distribuição de uma espécie e riqueza de espécies.
Esse capítulo segue parte da estrutura organizada por @lovelace-etal-2019, principalmente os capítulos 2 a 8, sendo adaptado para atender aos principais requisitos que julgamos necessários a estudos ecológicos. Entretanto, não foi possível cobrir todos os assuntos sobre geoprocessamento, sendo um assunto muito extenso, que requer a leitura de livros especializados na área como Wegmann et al. (2016), Wegmann et al. (2020) e Fletcher & Fortin (2018). Outros livros sobre a análise geoespacial no R podem ser consultados no capítulo 11 - Geospatial do Big Book of R.
Todas as operações serão realizadas através da linguagem R, utilizando principalmente os pacotes: tidyverse (Wickham et al. 2019) para o formato tidyverse, here (Müller 2020) para diretórios, sf [@pebesma-2018] para dados vetoriais, raster [@hijmans-2020] para dados raster, rgdal para formatos geoespaciais (Bivand et al. 2021), rnaturalearth (South 2017) e geobr (Pereira & Goncalves 2020) para baixar dados vetoriais, e ggplot2 (Wickham 2016), ggspatial (Dunnington 2020), tmap (Tennekes 2018), mapview (Appelhans et al. 2020), leaflet (Cheng et al. 2021), e viridis [@garnier-2018] para a composição de mapas, dentro outros.
Pré-requisitos
Dessa forma, garanta que esses pacotes listados a seguir estejam instalados e carregados.
# instalar pacotes
install.packages(c("tidyverse",
"here",
"sf",
"raster",
"rgdal",
"spData",
"rnaturalearth",
"geobr",
"ggplot2",
"ggspatial",
"tmap",
"tmaptools",
"grid",
"mapview",
"leaflet",
"viridis"),
dep = TRUE)# carregar pacotes
library(tidyverse)
library(here)
library(sf)
library(raster)
library(rgdal)
library(spData)
library(rnaturalearth)
library(geobr)
library(ggplot2)
library(ggspatial)
library(tmap)
library(tmaptools)
library(grid)
library(mapview)
library(leaflet)
library(viridis)IMPORTANTE: Se você estiver utilizando MacOS ou Linux, a instalação dos pacotes listados acima pode não funcionar. Esses sistemas operacionais possuem “requisitos específicos do sistema” que são geralmente descritos no
README.mddos pacotes no GitHub. Entretanto, há várias instruções específicas que podem ser encontradas on-line. Além deste link, há posts sobre Linux e MacOS que podema auxiliar.
Vetor
Dados vetoriais representam informações geográficas acuradas através de pontos, linhas e polígonos (Figura @ref(fig:fig-vetor-tipos)). Cada uma dessas geometrias são indicadas para representar feições e/ou eventos específicos, como veremos adiante.
Ilustração das geometrias de ponto, linha e polígono genéricos. Adaptado de: Lovelace, Nowosad & Muenchow (2019).
Pontos
Pontos são geometrias geralmente utilizadas para representar eventos pontuais unitários, como ocorrência de espécies, locais de coleta, pontos de GPS ou nascentes de rios. Esses dados são representados por um único vértice, ou seja, um par de coordenadas (longitude - X e latitude - Y), que são plotados na forma de cículos ou outro elemento que represente o evento em questão. Dessa forma, geralmente utilizamos dados tabulares com pelo menos duas colunas contendo essas coordenadas. Além disso, esses dados tabulares podem conter outras colunas com informações quantitativas ou qualitativas como número de espécies, temperatura, precipitação ou ainda categorias como tipo de habitat, que podemos representar nos pontos através de diferentes formatos, tamanhos ou cores desses pontos (Tabela @ref(tab:tab-vetor-pontos) e Figura @ref(fig:fig-vetor-pontos)).
| Id | Longitude | Latitude | Número de espécies | Temperatura | Precipitação | Habitat |
|---|---|---|---|---|---|---|
| 1 | 0 | 2 | 2 | 20 | 1000 | floresta |
| 2 | 1 | 5 | 3 | 22 | 1100 | pastagem |
| 3 | 2 | 3 | 3 | 28 | 1300 | floresta |
| 4 | 5 | 4 | 2 | 23 | 1200 | floresta |
| 5 | 5 | 1 | 5 | 25 | 1450 | pastagem |
Geometrias de pontos e suas identificações com a tabela de dados.
Linhas
Linhas representam geometrias lineares como estradas, rios, trajetos, divisões ou distâncias. Geralmente as linhas são criadas em softwares de Sistema de Informações Geográficas (SIG) como o QGIS, e depois importadas para o R. As linhas são representadas por no mínimo dois vértices conectados, i.e., dois pares de coordenadas, gerando uma geometria aberta, possuindo como característica o comprimento. Da mesma forma que os pontos, as linhas podem possuir informações tabulares, sendo quantitativas como comprimento dessa feição, ou ainda informações qualitativas como o nome de estradas ou vazão de rios, que podem ser utilizadas para alterar o formato, tamanho ou cor dessas linhas (Tabela @ref(tab:tab-vetor-linhas) e Figura @ref(fig:fig-vetor-linhas)).
| Id | Rodovias | Comprimento |
|---|---|---|
| 1 | rodovia_01 | 12 |
| 2 | rodovia_02 | 52 |
| 3 | rodovia_03 | 5 |
| 4 | rodovia_04 | 38 |
| 5 | rodovia_05 | 18 |
Geometrias de linhas e suas identificações com a tabela de dados.
Polígonos
Por fim, polígonos representam geometrias fechadas, como fragmentos de vegetação, lagos ou limites geográficos, sendo mais voltado para representar feições de um mapa de uso e cobertura da terra ou limites geográficos naturais, políticos, administrativos ou regulares. Os polígonos também são criados geralmente em softwares específicos de SIG e depois importados para o R, ou podemos usar funções para criar buffers ou malhas de quadrículas ou hexágonos. Os polígonos são representados por no mínimo três vértices conectados, sendo que o primeiro vértice possui coordenadas idênticas ao último, de modo que essa ligação gere uma feição fechada, com características como perímetro e/ou área. Da mesma forma que os pontos e linhas, colunas podem ser associadas aos polígonos para representar informações quantitativas como perímetro e área dessa polígono, ou ainda informações qualitativas como a classe de cobertura da terra ou o nome do limite geográfico, que podem ser utilizados para alterar formatos, tamanho ou cores desses polígonos (Tabela @ref(tab:tab-vetor-poligonos) e Figura @ref(fig:fig-vetor-poligonos)).
| id | uso | area_ha | perimeto_m |
|---|---|---|---|
| 1 | floresta | 50 | 700 |
| 2 | urbano | 22 | 300 |
| 3 | pastagem | 30 | 250 |
| 4 | agua | 25 | 400 |
| 5 | cerrado | 40 | 500 |
Geometrias de polígonos e suas identificações com a tabela de dados.
Além disso, geralmente utilizamos polígono regulares (buffers, quadrículas ou hexágonos) para resumir informações de biodiversidade ou de variáveis preditoras, que podem ser utilizadas como unidades amostrais em análises espaciais ou estatísticas, principalmente nas áreas de Ecologia Espacial, Ecologia da Paisagem, Biogeografia e Macroecologia (Tabela @ref(tab:tab-vetor-poligonos-regulares) e Figura @ref(fig:fig-vetor-poligonos-regulares)).
| id | numero_especies | temperatura | precipitacao |
|---|---|---|---|
| 1 | 2 | 20 | 1000 |
| 2 | 3 | 22 | 1100 |
| 3 | 3 | 28 | 1300 |
| 4 | 2 | 23 | 1200 |
| 5 | 5 | 25 | 1450 |
| … | … | … | … |
Polígonos regulares: buffers, quadrículas e hexágonos.
Tabela de atributos
Para os dados vetoriais é necessário ainda destacar um elemento fundamental: a tabela de atributos. A tabela de atributos é uma tabela que inclui dados geográficos e dados alfanuméricos. Os dados geográficos são representados por cada feição geolocalizada espacialmente (ponto, linha ou polígono), e os dados alfanuméricos são todos os demais dados associados a cada uma dessas feições, representado na forma de colunas (Figuras @ref(fig:fig-vetor-pontos), @ref(fig:fig-vetor-linhas), @ref(fig:fig-vetor-poligonos) e @ref(fig:fig-vetor-poligonos-regulares)).
Dessa forma, a tabela de atributos reúne informações sobre cada feição e pode ser utilizada para realizar de filtros ou agregações dos dados de cada feição. É nessa tabela que podemos ainda concatenar novas informações (colunas) de operações com as feições (linhas da tabela de atributos) como cálculo de comprimento, perímetro, área ou ainda outras operações com as colunas. Também podemos associar outros dados não espaciais aos dados da tabela de atributos com a junção por uma coluna identificadora.
sf: principal pacote no R para dados vetoriais
Atualmente o principal pacote para trabalhar com dados vetoriais no R é o sf, que implementou o Simple Feature no R (Pebesma 2018). Entretanto, outro pacote pode ser tão versátil quanto o sf, no caso o terra, ainda em desenvolvimento.
O pacote sf facilitou muito a forma de trabalhar com vetores no R, sendo que as principais vantagens desse pacote são (Lovelace, Nowosad & Muenchow 2019):
- rápida importação e exportação de dados
- aprimorado desempenho de plotagem
- objetos sf podem ser tratados como dataframes na maioria das operações
- funções sf podem ser combinadas usando o operador
%>%e funcionam no formato tidyverse - funções sf são consistentes e intuitivas (quase sempre começam com prefixo
st_)
Os tipos de geometrias apresentadas são representadas por diferentes classes: POINT, LINESTRING e POLYGON para apenas uma feição de cada tipo de geometria; MULTIPOINT, MULTILINESTRING e MULTIPOLYGON para várias feições de cada tipo de geometria e; GEOMETRYCOLLECTION para várias feições e tipos de geometrias (Figura @ref(fig:fig-vetor-sf-classes)).
Tipos de classes suportadas pelo pacote sf. Fonte: Lovelace, Nowosad & Muenchow (2019).
O pacote sf define um sistema de três classes hierárquicas (Tabela @ref(tab:tab-vetor-sf-estruturas)):
- Classe sfg - uma geometria única
- Classe sfc - uma coluna de geometria, que é um conjunto de geometrias sfg e informações Sistema de Referência de Coordenadas (do inglês Coordinate Reference Systems- CRS)
- Classe sf - uma camada, que é uma coluna de geometria sfc dentro de um dataframe com atributos não espaciais
| Classes | Hierarquia | Informação |
|---|---|---|
sfg |
Geometria | Tipo e coordenadas |
sfc |
Coluna de geometria | Conjundo de sfg + CRS |
| sf | Camada | sfc + atributos |
Ao olharmos as informações de um objeto da classe sf, podemos notar diversas informações que descrevem o mesmo:
- resumo do vetor: indica o número de feições (linhas) e campos (colunas)
- tipo da geometria: umas das sete classes listadas anteriormente (Figura @ref(fig:fig-vetor-sf-classes))
- dimensão: número de dimensões, geralmente duas (XY)
- bbox (bordas): coordenadas mínimas e máximas da longitude e latitude
- informação do CRS:
epsgouproj4stringindicando o CRS - tibble: tabela de atributos, com destaque para a coluna
geomougeometryque representa cada feição ou geometria
data(world)
world## Simple feature collection with 177 features and 10 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## geographic CRS: WGS 84
## # A tibble: 177 x 11
## iso_a2 name_long continent region_un subregion type area_km2 pop lifeExp
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 FJ Fiji Oceania Oceania Melanesia Sove… 1.93e4 8.86e5 70.0
## 2 TZ Tanzania Africa Africa Eastern … Sove… 9.33e5 5.22e7 64.2
## 3 EH Western … Africa Africa Northern… Inde… 9.63e4 NA NA
## 4 CA Canada North Am… Americas Northern… Sove… 1.00e7 3.55e7 82.0
## 5 US United S… North Am… Americas Northern… Coun… 9.51e6 3.19e8 78.8
## 6 KZ Kazakhst… Asia Asia Central … Sove… 2.73e6 1.73e7 71.6
## 7 UZ Uzbekist… Asia Asia Central … Sove… 4.61e5 3.08e7 71.0
## 8 PG Papua Ne… Oceania Oceania Melanesia Sove… 4.65e5 7.76e6 65.2
## 9 ID Indonesia Asia Asia South-Ea… Sove… 1.82e6 2.55e8 68.9
## 10 AR Argentina South Am… Americas South Am… Sove… 2.78e6 4.30e7 76.3
## # … with 167 more rows, and 2 more variables: gdpPercap <dbl>,
## # geom <MULTIPOLYGON [°]>Podemos fazer um mapa simples utilizando a função plot() desse objeto. Para facilitar, escolheremos apenas a primeira coluna [1] (Figura @ref(fig:fig-vetor-mundo)).
IMPORTANTE: faremos mapas mais elaborados na seção xx desse capítulo.
plot(world[1], col = viridis::viridis(100), main = "Mapa do mundo")## Warning in plot.sf(world[1], col = viridis::viridis(100), main = "Mapa do
## mundo"): col is not of length 1 or nrow(x): colors will be recycled; use pal to
## specify a color palette
Mapa vetorial do mundo.
Raster
Os dados no formato raster consistem em uma matriz (com linhas e colunas) representando células igualmente espaçadas (pixels; Figura @ref(fig:fig-raster)). Esse formato de dado torna a álgebra e o processamento de mapas muito mais eficiente e rápido do que o processamento de dados vetoriais. As células dos dados raster possuem duas informações: 1. identificação das células (IDs das células) para especificar sua posição na matriz (Figura @ref(fig:fig-raster) A) e; 2. valores das células (Figura @ref(fig:fig-raster) B), que geralmente são coloridos para facilitar a interpretação da variação dos valores no espaço (Figura @ref(fig:fig-raster) C). Além disso, valores ausentes ou não amostrados são representados por NA, ou seja, not available (Figura @ref(fig:fig-raster) B e C).
##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combine
Raster: (A) IDs das células, (B) valores das células, (C) células coloridas. Adaptado de: Lovelace, Nowosad & Muenchow (2019).
Tipos de raster
A célula ou pixel de um raster pode conter apenas um único valor, que pode ser contínuo ou categórico (Figura @ref(fig:fig-raster-cont-cat)). O formato raster geralmente representa fenômenos contínuos, como elevação, precipitação, temperatura, ou dados espectrais de imagens de satélite, mas também pode representar categorias como tipos de florestas ou cobertura da terra (Figura @ref(fig:fig-raster-cont-cat)).
Raster: (A) mapa contínuo, (B) mapa categórico. Adaptado de: Lovelace, Nowosad & Muenchow (2019).
raster: principal pacote no R para dados raster
Atualmente, o principal pacote para trabalhar com dados raster é o raster (Hijmans 2020), apesar de existir outros dois em desenvolvimento e já sendo aplicados, como o terra e o stars. O pacote raster fornece uma ampla gama de funções para criar, importar, exportar, manipular e processar dados raster no R. O objeto raster pode assumir três classes no R: RasterLayer, RasterStack e RasterBrick.
A classe RasterLayer representa apenas uma camada raster. Para criar um raster no R podemos utilizar a função raster::raster(). Observando essa classe, podemos notar as seguintes informações:
- class: classe raster do objeto raster
- dimensions: número de linhas, colunas e células
- resolution: largura e altura da célula
- extent: coordenadas mínimas e máximas da longitude e latitude
- crs: Sistema de Referência de Coordenadas (CRS)
- source: fonte dos dados (memória ou disco)
- names: nome das camadas
- values: valores máximos e mínimos das células
raster_layer <- raster::raster(nrows = 5, ncols = 5,
res = .5,
xmn = -61.5, xmx = -59, ymn = -14.5, ymx = -12,
vals = sample(1:25, 25, rep = TRUE))
raster_layer## class : RasterLayer
## dimensions : 5, 5, 25 (nrow, ncol, ncell)
## resolution : 0.5, 0.5 (x, y)
## extent : -61.5, -59, -14.5, -12 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 1, 24 (min, max)Um mapa simples do objeto raster pode ser obtido utilizando a função plot(), do próprio pacote raster (Figura @ref(fig:fig-raster-layer)).
plot(raster_layer, col = viridis::viridis(n = 25))
Mapa simples de um RasterLayer.
Além da classe RasterLayer, há mais duas classes que trabalham com múltiplas camadas: RasterBrick e RasterStack. Elas diferem em relação ao número de formatos de arquivo suportados, tipo de representação interna e velocidade de processamento.
A classe RasterBrick geralmente corresponde a importação de um único arquivo de imagem de satélite multiespectral (multicamadas) ou a um único objeto com várias camadas na memória. A função raster::brick() cria um objeto RasterBrick.
raster_layer1 <- raster_layer
raster_layer2 <- raster_layer * raster_layer
raster_layer3 <- sqrt(raster_layer)
raster_layer4 <- log10(raster_layer)
raster_brick <- raster::brick(raster_layer1, raster_layer2, raster_layer3, raster_layer4)
raster_brick## class : RasterBrick
## dimensions : 5, 5, 25, 4 (nrow, ncol, ncell, nlayers)
## resolution : 0.5, 0.5 (x, y)
## extent : -61.5, -59, -14.5, -12 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer.1, layer.2, layer.3, layer.4
## min values : 1, 1, 1, 0
## max values : 24.000000, 576.000000, 4.898979, 1.380211Ao utilizarmos a função plot() do pacote raster, podemos visualizar todos os raster contidos no objeto RasterBrick (Figura @ref(fig:fig-raster-brick)).
plot(raster_brick, col = viridis::viridis(n = 25))
Mapas simples de um raster RasterBrick.
Já a classe RasterStack permite conectar vários objetos raster armazenados em arquivos diferentes ou vários objetos na memória. Um RasterStack é uma lista de objetos RasterLayer com a mesma extensão, resolução e CRS. Uma maneira de criá-lo é com a junção de vários objetos espaciais já existentes no ambiente do R ou listar vários arquivos raster em um diretório armazenado no disco. A função raster::stack() cria um objeto RasterStack.
Outra diferença é que o tempo de processamento para objetos RasterBrick geralmente é menor do que para objetos RasterStack. A decisão sobre qual classe Raster deve ser usada depende principalmente do caráter dos dados de entrada.
raster_layer1 <- raster_layer
raster_layer2 <- raster_layer * raster_layer
raster_layer3 <- sqrt(raster_layer)
raster_layer4 <- log10(raster_layer)
raster_stack <- raster::stack(raster_layer1, raster_layer2, raster_layer3, raster_layer4)
raster_stack## class : RasterStack
## dimensions : 5, 5, 25, 4 (nrow, ncol, ncell, nlayers)
## resolution : 0.5, 0.5 (x, y)
## extent : -61.5, -59, -14.5, -12 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : layer.1, layer.2, layer.3, layer.4
## min values : 1, 1, 1, 0
## max values : 24.000000, 576.000000, 4.898979, 1.380211Da mesma forma, ao utilizar a função plot() do pacote raster, podemos visualizar todos os raster contidos no objeto RasterStack (Figura @ref(fig:fig-raster-stack)).
plot(raster_stack, col = viridis::viridis(n = 25))Mapas simples de um raster RasterStack.
Sistema de Referência de Coordenadas e Unidades
Os dados geográficos (vetor e raster) possuem ainda um outro componente fundamental, que é o Sistema de Referência de Coordenadas, ou do inglês Coordinate Reference System (CRS). Esse componente define como os elementos espaciais (vetor e raster) representam uma feição na superfície da Terra. Esse componente é composto por dois principais conceitos: primeiro, que tipo de unidades estão sendo utilizadas para a representação geográfica, podendo assumir dois tipos - ângulos ou metros, que definem o Sistema de Coordenadas Geográficas e o Sistema de Coordenadas Projetadas, respectivamente.
O segundo componente é o datum, que é a relação do sistema de coordenadas com a superfície da Terra. Esse último componente faz parte de uma área da Cartografia denominada Geodésia que estuda a forma e dimensões da Terra, campo gravitacional, e a localização de pontos fixos e sistemas de coordenadas. O livro de Lapaine et al. (2017) é um excelente material para se aprofundar nesse assunto.
Sistema de Coordenadas Geográficas
O Sistema de Referência de Coordenadas Geográficas utiliza ângulos para representar feições na superfície da Terra através de dois valores: longitude e latitude. A longitude localiza-se na direção Leste-Oeste e a latitude localiza-se na direção Norte-Sul. Nesse sistema, a superfície da Terra geralmente é representada por uma superfície elipsoidal, pois a Terra é ligeiramente achatada nos polos.
Sistema de Coordenadas Projetadas
O Sistema de Referência de Coordenadas Projetadas utiliza um Sistema Cartesiano de Coordenadas em uma superfície plana. Dessa forma, à partir de uma origem, traçam-se eixos x e y, e uma unidade linear como o metro é utilizada. Todos as projeções feitas de sistemas geográficos convertem uma superfície tridimensional em uma superfície plana bidimensional. Sendo assim, essa conversão trás consigo algum tipo de distorção em relação à porção real, podendo ser distorções em: 1. formas locais, 2. áreas, 3. distâncias , 4. flexão ou curvatura , 5. assimetria, e 6. lacunas de continuidade. Dessa forma, um sistema de coordenadas projetadas pode preservar somente uma ou duas dessas propriedades.
Existem três grandes grupos de projeções: cilíndricos, cônicos e planares. Na projeção cilíndrica, a superfície da Terra é mapeada em um cilindro, sendo também criada tocando a superfície da Terra ao longo de uma ou duas linhas de tangência, sendo utilizada com mais frequência para mapear todo o globo, tendo como exemplo mais conhecido a Projeção Universal Transversa de Mercator (UTM). Na projeção cônica, a superfície da Terra é projetada em um cone ao longo de uma linha ou duas linhas de tangência, de modo que as distorções são minimizadas ao longo das linhas e aumentam com a distância das mesmas, sendo portanto, mais adequada para mapear áreas de latitudes médias, tendo como exemplo mais conhecido a Projeção Cônica Equivalente de Albers e a Projeção Cônica Conforme de Lambert. E na projeção plana, também denominada Projeção azimutal, o mapeamento toca o globo em um ponto ou ao longo de uma linha de tangência, sendo normalmente usado no mapeamento de regiões polares, sendo a mais comum a Projeção Azimutal Equidistante, a mesma utilizada na bandeira da ONU.
Datum
Como dito anteriormente, o datum é a relação do sistema de coordenadas com a superfície da Terra. Ele representa o ponto de intersecção do elipsoide de referência com a supercie da Terra (geoide, forma verdadeira da Terra), compensando as diferenças do campo gravitacional da Terra. Existem dois tipos de datum - local e geocêntrico. Em um datum local, como o SAD69 - South American Datum 1969, o elipsoide de referência é deslocado para se alinhar com a superfície em um determinado local, por exemplo, na América do Sul. Já em um datum geocêntrico, como WGS84 - World Geodetic System 1984, o centro do elipsoide é o centro de gravidade da Terra e a precisão das projeções não é otimizada para um local específico do globo.
No Brasil, desde 2015, o Instituto Brasileiro de Geografia e Estatística (IBGE) adotou utilizar o datum SIRGAS2000 - Sistema de Referencia Geocéntrico para las Américas 2000 para todos os mapeamentos realizados no Brasil, um esforço conjunto para adotar o mesmo datum em toda a América. Mais sobre esse datum pode ser lido aqui: SIRGAS2000.
Sistema de Referência de Coordenadas (CRS) no R
No R, há duas formas principais de representar um Sistema de Referência de Coordenadas: 1. código epsg e 2. proj4string. O código EPSG (European Petroleum Survey Group) é uma sequência de números curta, referindo-se apenas a um CRS. O site epsg.io permite consultar diversas informações sobre um código, como procurar por um código, representação de mapas e fazer transformações de CRS.
Já proj4string permite mais flexibilidade para especificar diferentes parâmetros, como o tipo de projeção, datum e elipsóide. Dessa forma, é possível especificar muitas projeções, ou mesmo modificar as projeções existentes, tornando a representação proj4string mais complexa e flexível.
Além disso, ainda é possível consultar uma extensa lista de CRSs no site spatialreference.org, que fornece descrições em diversis formatos, baseados em GDAL e Proj.4. Essa abordagem permite consultar uma URL que pode produzir uma referência espacial em um formato que seu software SIG ou o R pode utilizar como referência.
Os pacotes espaciais no R suportam uma ampla variedade de CRSs e usam a biblioteca PROJ. A função rgdal::make_EPSG() retorna um data frame das projeções disponíveis, com informações dos códigos epsg e proj4string numa mesma tabela, facilitando a busca e uso de CRSs (Tabela @ref(tab:tab-epsg)).
crs_data <- rgdal::make_EPSG()
head(crs_data)| code | note | prj4 | prj_method |
|---|---|---|---|
| 3819 | HD1909 | +proj=longlat +ellps=bessel +no_defs +type=crs | (null) |
| 3821 | TWD67 | +proj=longlat +ellps=aust_SA +no_defs +type=crs | (null) |
| 3822 | TWD97 | +proj=geocent +ellps=GRS80 +units=m +no_defs +type=crs | (null) |
| 3823 | TWD97 | +proj=longlat +ellps=GRS80 +no_defs +type=crs | (null) |
| 3824 | TWD97 | +proj=longlat +ellps=GRS80 +no_defs +type=crs | (null) |
| 3887 | IGRS | +proj=geocent +ellps=GRS80 +units=m +no_defs +type=crs | (null) |
Principais fontes de dados geográficos
Existe diversas fontes de dados geográficos em diferentes bases de dados disponíveis gratuitamente. Geralmente essas bases de dados são disponibilizadas separadamente em apenas dados vetoriais e dados raster. Para dados vetoriais, grande parte dos dados disponibilizados são utilizados em mapas como limites políticos, limites de biomas ou distribuição de espécies para polígonos; estradas e rios para dados lineares, ou ainda pontos de ocorrência de espécies ou comunidades, ou medidas tomadas em campo sobre condições naturais como clima ou relevo, como pontos. Entretanto, é sempre recomendado o uso de bases oficiais, principalmente em relação a dados vetoriais de limites políticos. Para tanto, é fundamental buscar as bases oficiais de cada país, entrentanto, há bases que podem ser utilizadas globalmente, como veremos.
Sobre as bases de dados raster, há uma infinidade de dados para diferentes objetivos, mas grande parte deles são relativos à condições ambientais, representando uma variável de interesse de forma contínua no espaço, como temperatura, precipitação, elevação, etc.
Há uma compilação de dados geográficos vetoriais e raster feita por Marcus Vinícius Alves de Carvalho e Angelica Carvalho Di Maio, chamada GeoLISTA. Entretanto, como as bases de dados tendem à ser muito dinâmicas, é possível que muitas bases tenham surgido e desaparecido desde a listagem realizada.
Além das bases de dados, há pacotes específicos no R que fazem o download de dados vetoriais e rasters, facilitando a aquisição e reprodutibilidade. Para conferir uma listagem completa de pacotes para diversas análises espaciais, veja CRAN Task View: Analysis of Spatial Data.
Vetor
Dentre as bases vetoriais, destacamos as seguintes na Tabela @ref(tab:tab-vetor-bases):
| Bases de dados | Descrição |
|---|---|
| IBGE | Limites territoriais e censitários do Brasil |
| FBDS | Uso da terra, APP e hidrografia - Mata Atlântica e Cerrado |
| GeoBank | Dados geológicos do Brasil |
| Pastagem.org | Dados de pastagens e gado para o Brasil |
| CanaSat | Dados de cana-de-açúcar para o Brasil |
| CSR Maps | Diversos dados vetoriais e raster para o Brasil |
| Ecoregions | Dados de biorregiões e biomas do mundo |
| UN Biodiversity Lab | Diversas bases de dados para o mundo |
| Biodiversity Hotspots | Dados dos limites dos Hotspots de Biodiversidade |
| IUCN Red List of Threatened Species | Dados dos limites das distribuições das espécies para o mundo |
| Map of Life (MOL) | Dados da distribuição de espécies e outros dados para o mundo |
| Key Biodiversity Areas | Dados dos limites das Key Biodiversity Areas |
| HydroSHEDS | Informações hidrológicas do mundo |
| Global Roads Inventory Project (GRIP) | Dados de estradas do mundo todo |
| Database of Global Administrative Areas (GADM) | Limites de áreas administrativas do mundo |
| Natural Earth | Diversos limites para o mundo |
| Protected Planet | Limites de áreas protegidas para o mundo |
| Global Biological Information Facility (GBIF) | Dados de ocorrências de espécies para o mundo |
| Species Link | Dados de ocorrências de espécies para o Brasil |
| Global Invasive Species Information Network (GISIN) | Dados de ocorrências de espécies invasoras para o Mundo |
Raster
Dentre as bases raster, destacamos as seguintes na Tabela @ref(tab:tab-raster-bases):
| Bases de dados | Descrição |
|---|---|
| MapBiomas | Uso e cobertura da terra para o Brasil, Panamazonia Legal e Chaco, de 1985 a 2019 |
| Bahlu | Distribuições históricas de terras agrícolas e pastagens para todo o Brasil de 1940 a 2012 |
| USGS | Dados de diversos satélites livres para o mundo |
| SRTM | Dados de elevação para o mundo |
| Geoservice Maps | Dados de elevação e florestas para o mundo |
| Global Forest Watch | Dados de florestas para o mundo |
| GlobCover | Dados de uso e cobertura da terra para todo o planeta |
| Landcover | Dados de uso e cobertura da terra para todo o planeta |
| Global Human Footprint | Dados de pegada ecológica para o mundo |
| GHSL - Global Human Settlement Layer | Dados e ferramentas abertos e gratuitos para avaliar a presença humana no planeta |
| Land-Use Harmonization (LUH2) | Dados atuais e previsões de uso da terra |
| ESA Climate Change Initiative | Arquivos globais de observação da Terra nos últimos 30 anos da Agência Espacial Europeia (ESA) |
| WorldClim | Dados climáticos para o mundo |
| CHELSA | Dados climáticos para o mundo |
| EarthEnv | Dados de cobertura da terra, nuvens, relevo e hidrografia |
| SoilGrids | Dados de solo para o mundo |
| Global Wetlands | Dados de áreas úmidas para o mundo |
| Global Surface Water Explorer | Dados de águas superficiais para o mundo |
| MARSPEC | Dados de condições do oceano para o mundo |
| Bio-ORACLE | Dados de condições do oceano para o mundo |
Pacotes do R
Dentre os pacotes no R para download de dados geográficos, destacamos os seguintes na Tabela @ref(tab:tab-packages-bases):
| Pacotes | Descrição |
|---|---|
| geobr | Carrega Shapefiles de Conjuntos de Dados Espaciais Oficiais do Brasil |
| rnaturalearth | Dados do mapa mundial da Natural Earth |
| rworldmap | Mapeando Dados Globais |
| spData | Conjuntos de dados para análise espacial |
| OpenStreetMap | Acesso para abrir imagens raster de mapas de ruas |
| osmdata | Baixe e importe dados do OpenStreetMap |
| geonames | Interface para o serviço da Web de consulta espacial ‘Geonames’ |
| rgbif | Interface para o Global ‘Biodiversity’ Information Facility API |
| maptools | Ferramentas para lidar com objetos espaciais |
| marmap | Importar, traçar e analisar dados batimétricos e topográficos |
| oce | Fonte e processamento de dados oceanográficos |
| envirem | Geração de Variáveis ENVIREM |
| sdmpredictors | Conjuntos de dados preditor de modelagem de distribuição de espécies |
| metScanR | Encontre, Mapeie e Colete Dados e Metadados Ambientais |
| ClimDown | Biblioteca de redução de escala do clima para a produção diária do modelo climático |
| rWBclimate | Acessa dados climáticos do Banco Mundial |
| rnoaa | Dados meteorológicos ‘NOAA’ de R |
| RNCEP | Obtenha, organize e visualize dados meteorológicos NCEP |
| smapr | Aquisição e processamento de dados ativos-passivos (SMAP) de umidade do solo da NASA |
Importar e exportar dados geográficos
Agora que sabemos o que são dados geográficos e em quais bases de dados podemos buscar e baixar esses dados, veremos seus principais formatos e como importá-los e exportá-los do R.
Principais formatos de arquivos geográficos
Há diversos formatos de arquivos geográficos, alguns específicos para dados vetoriais e raster, e outros no formato de banco de dados geoespaciais, como PostGIS, que podem armazenar ambos os formatos.
Entretanto, todos os formatos para serem importados para o R usam do GDAL (Geospatial Data Abstraction Library), uma interface unificada para leitura e gravação de diversos formatos de arquivos geográficos, sendo utilizado também por uma séria de softwares de GIS como QGIS, GRASS GIS e ArcGIS.
Dentre esses formatos, destacamos os seguintes na Tabela @ref(tab:tab-formatos).
| Nome | Extenção | Descrição | Tipo | Modelo |
|---|---|---|---|---|
| ESRI Shapefile | .shp (arquivo principal) | Formato popular que consiste em pelo menos quatro arquivos: .shp (feição), .dbf (tabela de atributos), .shx (ligação entre .shp e .dbf) e .prj (projeção) | Vetor | Parcialmente aberto |
| GeoJSON | .geojson | Estende o formato de troca JSON incluindo um subconjunto da representação de recurso simples | Vetor | Aberto |
| KML | .kml | Formato baseado em XML para visualização espacial, desenvolvido para uso com o Google Earth. O arquivo KML compactado forma o formato KMZ | Vetor | Aberto |
| GPX | .gpx | Esquema XML criado para troca de dados de GPS | Vetor | Aberto |
| GeoTIFF | .tif/.tiff | Formato raster popular. Um arquivo TIFF contendo metadados espaciais adicionais. | Raster | Aberto |
| Arc ASCII | .asc | Formato de texto em que as primeiras seis linhas representam o cabeçalho raster, seguido pelos valores das células raster organizadas em linhas e colunas | Raster | Aberto |
| NetCDF | .nc | NetCDF (Network Common Data Form) é um conjunto de bibliotecas de software e formatos de dados independentes para criação | Raster | Aberto |
| BIL | .bil/.hdr | BIL (Banda intercalada por linha) são métodos comuns de organização para imagens multibanda, geralmente acompanhados por um arquivo .hdr, descrevendo atributos específicos da imagem | Raster | Aberto |
| R-raster | .gri/ .grd | Formato raster nativo do raster do pacote R | Raster | Aberto |
| SQLite/SpatiaLite | .sqlite | Banco de dados relacional autônomo | Vetor e raster | Aberto |
| ESRI FileGDB | .gdb | Objetos espaciais e não espaciais criados pelo ArcGIS. Permite: várias classes de recursos; topologia | Vetor e raster | Proprietário |
| GeoPackage | .gpkg | Contêiner de banco de dados leve baseado em SQLite permitindo uma troca fácil e independente de plataforma de geodados | Vetor e raster | Aberto |
O formato mais comum para arquivos vetoriais é o ESRI Shapefile, para arquivos raster é o GeoTIFF, e para dados climáticos em múltiplas camadas, geralmente há a disponibilização de dados no formato NetCDF. Entretanto, recentemente tivemos o surgimento do GeoPackage, que possui diversas vantagens em relação aos formatos anteriores, podendo armazanar em apenas um arquivo, dados no formato vetorial, raster e também dados não-espaciais, além de possuir uma grande integração com diversos softwares e bancos de dados.
Importar dados
As principais funções para importar dados no R são: 1) para vetores a função sf::st_read(), e 2) para raster a função raster::raster(), e suas variações raster::brick() e raster::stack() para múltiplas camadas. Essas funções atribuem objetos ao seu espaço de trabalho, armazenando-os na memória RAM disponível em seu hardware, sendo essa a maior limitação para trabalhar com dados geográficos no R. Por exemplo, se um arquivo raster possui mais de 8 Gb de tamamho, e seu computador possui extamente 8 Gb de RAM, é muito provável que ele não seja importado ou mesmo criado como um objeto dentro do ambiente R. Existem soluções para esses problemas, mas não as abordaremos nesse capítulo.
Vetor
Como vimos, os arquivos vetoriais são disponibilizados em diversos formatos. Para sabermos se um determinado formato pode ser importado ou exportado utilizando o pacote sf, podemos utilizar a função sf::st_drivers(). Uma amostra desses formatos é apresentado na Tabela @ref(tab:tab-vetor-formatos):
head(sf::st_drivers())| name | long_name | write | copy | is_raster | is_vector | vsi |
|---|---|---|---|---|---|---|
| ESRIC | Esri Compact Cache | FALSE | FALSE | TRUE | TRUE | TRUE |
| FITS | Flexible Image Transport System | TRUE | FALSE | TRUE | TRUE | FALSE |
| PCIDSK | PCIDSK Database File | TRUE | FALSE | TRUE | TRUE | TRUE |
| netCDF | Network Common Data Format | TRUE | TRUE | TRUE | TRUE | TRUE |
| PDS4 | NASA Planetary Data System 4 | TRUE | TRUE | TRUE | TRUE | TRUE |
| VICAR | MIPL VICAR file | TRUE | TRUE | TRUE | TRUE | TRUE |
Importar dados vetoriais existentes
Para importar vetores existentes para o R, iremos utilizar a função sf::st_read(). A estrutura é semelhante para todos os formatos descritos na Tabela @ref(tab:tab-vetor-formatos), de modo que sempre preencheremos o argumento dsn (data source name) com o nome do arquivo a ser importado. Entretanto, para banco de dados, como GeoPackage, pode ser necessário especificar a camada que se tem interesse com um segundo argumento chamado layer, com o nome da camada.
Para quase todas as operações vetoriais nesse capítulo, usaremos os dados disponíveis para o município de Rio Claro/SP. Primeiramente, iremos baixar esses dados da FBDS (Fundação Brasileira para o Desenvolvimento Sustentável), através desse repositório de dados. Em 2013, a FBDS deu início ao Projeto de Mapeamento em Alta Resoluçăo dos Biomas Brasileiros, mapeando a cobertura da terra, hidrografia (nascentes, rios e lagos) e áreas de preservaçăo permanente (APPs). O mapeamento foi concluído para os municípios dos biomas Mata Atlântica e Cerrado. Para fazer o download dos arquivos de interesse, utilizaremos o R, através da função download.file().
Primeiramente, iremos criar um diretório com a função create.dir(), usando a função here::here() para indicar o repositório (ver seção xx).
# criar um diretorio
dir.create(here::here("dados"))
dir.create(here::here("dados", "vetor"))Em seguida, vamos fazer o download de pontos de nascentes, linhas de hidrografia e polígonos de cobertura da terra para o município de Rio Claro/SP.
# aumentar o tempo de download
options(timeout = 600)
# download
for(i in c(".dbf", ".prj", ".shp", ".shx")){
# pontos de nascentes
download.file(
url = paste0("http://geo.fbds.org.br/SP/RIO_CLARO/HIDROGRAFIA/SP_3543907_NASCENTES", i),
destfile = here::here("dados", "vetor", paste0("SP_3543907_NASCENTES", i)), mode = "wb")
# linhas de hidrografia
download.file(
url = paste0("http://geo.fbds.org.br/SP/RIO_CLARO/HIDROGRAFIA/SP_3543907_RIOS_SIMPLES", i),
destfile = here::here("dados", "vetor", paste0("SP_3543907_RIOS_SIMPLES", i)), mode = "wb")
# poligonos de cobertura da terra
download.file(
url = paste0("http://geo.fbds.org.br/SP/RIO_CLARO/USO/SP_3543907_USO", i),
destfile = here::here("dados", "vetor", paste0("SP_3543907_USO", i)), mode = "wb")
}Agora podemos importar esses dados para o R. Primeiro vamos importar as nascentes (Figura @ref(fig:fig-vetor-nascentes)).
# importar pontos
rc_nas <- sf::st_read(here::here("dados", "vetor", "SP_3543907_NASCENTES.shp"), quiet = TRUE)plot(rc_nas[1], pch = 20, col = "blue", main = NA, axes = TRUE, graticule = TRUE)
Mapa de nascentes de Rio Claro/SP.
Agora vamos importar a hidrografia (Figura @ref(fig:fig-vetor-hidrografia)).
# importar hidrografia
rc_hid <- sf::st_read(here::here("dados", "vetor", "SP_3543907_RIOS_SIMPLES.shp"), quiet = TRUE)plot(rc_hid[1], col = "steelblue", main = NA, axes = TRUE, graticule = TRUE)
Mapa da hidrografia de Rio Claro/SP.
E por fim, vamos importar a cobertura da terra (Figura @ref(fig:fig-vetor-cobertura)).
# importar cobertura da terra
rc_cob <- sf::st_read(here::here("dados", "vetor", "SP_3543907_USO.shp"), quiet = TRUE)plot(rc_cob[5], col = c("blue", "orange", "gray30", "forestgreen", "green"), main = NA, axes = TRUE, graticule = TRUE)
Mapa de cobertura da terra de Rio Claro/SP.
Importar utilizando pacotes
Além de dados existentes, podemos importar dados vetoriais de pacotes, como listado anteriormente na Tabela @ref(tab:tab-packages-bases). Para o Brasil, o pacote mais interessante trata-se do geobr, do Instituto de Pesquisa Econômica Aplicada (IPEA), que possui dados oficiais do Instituto Brasileiro de Geografia e Estatística (IBGE)).
É possível listar todos os dados disponíveis no pacote através da função geobr::list_geobr(). Na Tabale @ref(tab:tab-vetor-dados-geobr) é possível ver alguns desses dados.
# listar todos os dados do geobr
geobr::list_geobr()| function | geography | years | source |
|---|---|---|---|
read_country |
Country | 1872, 1900, 1911, 1920, 1933, 1940, 1950, 1960, 1970, 1980, 1991, 2000, 2001, 2010, 2013, 2014, 2015, 2016, 2017, 2018, 2019 | IBGE |
read_region |
Region | 2000, 2001, 2010, 2013, 2014, 2015, 2016, 2017, 2018, 2019 | IBGE |
read_state |
States | 1872, 1900, 1911, 1920, 1933, 1940, 1950, 1960, 1970, 1980, 1991, 2000, 2001, 2010, 2013, 2014, 2015, 2016, 2017, 2018, 2019 | IBGE |
read_meso_region |
Meso region | 2000, 2001, 2010, 2013, 2014, 2015, 2016, 2017, 2018, 2019 | IBGE |
read_micro_region |
Micro region | 2000, 2001, 2010, 2013, 2014, 2015, 2016, 2017, 2018, 2019 | IBGE |
read_intermediate_region |
Intermediate region | 2017, 2019 | IBGE |
Como exemplo, vamos fazer o download o limite do município de Rio Claro/SP, utilizando o código do município (3543907) (Figura @ref(fig:fig-vetor-rio-claro)).
# rio claro
rc_2019 <- geobr::read_municipality(code_muni = 3543907, year = 2019, showProgress = FALSE)## Using year 2019plot(rc_2019[1], col = "gray", main = NA, axes = TRUE, graticule = TRUE)
Limite do município de Rio Claro/SP.
Já para o mundo, o pacote mais interessante trata-se do rnaturalearth, que faz o download de dados do Natural Earth. Vamos fazer o download do limite do Brasil (Figura @ref(fig:fig-vetor-brasil)).
# brasil
br <- rnaturalearth::ne_countries(scale = "large", country = "Brazil", returnclass = "sf")## Warning in fun(libname, pkgname): rgeos: versions of GEOS runtime 3.9.0-CAPI-1.16.2
## and GEOS at installation 3.8.0-CAPI-1.13.1differplot(br[1], col = "gray", main = NA, axes = TRUE, graticule = TRUE)
Limite do Brasil.
Criar um objeto espacial de uma tabela de coordenadas
É muito comum em coletas de campo ou fontes de dados, ter coordendas de locais de estudo ou de ocorrências de espécies organizadas em tabelas. Essas tabelas devem possuir duas colunas: longitude e latitude, ou X e Y para dados UTM, por exemplo. Ao importá-las para o R, o formato que assumem pode ser de uma das classes: matrix, dataframe ou tibble, ou seja, ainda não são da classe vetorial sf. Nesta seção iremos ver como fazer essa conversão.
Para tanto, vamos usar os dados de comunidades de anfíbios da Mata Atlântica (Atlantic Amphibians, Vancine et al. 2018). Iremos fazer o download diretamente do site da fonte dos dados. Antes vamos criar um diretório.
# criar um diretorio
dir.create(here::here("dados", "tabelas"))Em seguida, vamos fazer o download de um arquivo .zip e vamos extrair usando a função unzip() nesse mesmo diretório.
# download
download.file(url = "https://esajournals.onlinelibrary.wiley.com/action/downloadSupplement?doi=10.1002%2Fecy.2392&file=ecy2392-sup-0001-DataS1.zip",
destfile = here::here("dados", "tabelas", "atlantic_amphibians.zip"), mode = "wb")
# unzip
unzip(zipfile = here::here("dados", "tabelas", "atlantic_amphibians.zip"),
exdir = here::here("dados", "tabelas"))Agora podemos importar a tabela de dados com a função `readr::read_csv()``.
# importar tabela de locais
aa_locais <- readr::read_csv(here::here("dados", "tabelas", "ATLANTIC_AMPHIBIANS_sites.csv"))##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_character(),
## reference_number = col_double(),
## species_number = col_double(),
## month_start = col_double(),
## year_start = col_double(),
## month_finish = col_double(),
## year_finish = col_double(),
## effort_months = col_double(),
## latitude = col_double(),
## longitude = col_double(),
## altitude = col_double(),
## temperature = col_double(),
## precipitation = col_double()
## )
## ℹ Use `spec()` for the full column specifications.aa_locais## # A tibble: 1,163 x 25
## id reference_number species_number record sampled_habitat active_methods
## <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 amp1001 1001 19 ab fo,ll as
## 2 amp1002 1002 16 co fo,la,ll as
## 3 amp1003 1002 14 co fo,la,ll as
## 4 amp1004 1002 13 co fo,la,ll as
## 5 amp1005 1003 30 co fo,ll,br as
## 6 amp1006 1004 42 co tp,pp,la,ll,is <NA>
## 7 amp1007 1005 23 co sp as
## 8 amp1008 1005 19 co sp,la,sw as,sb,tr
## 9 amp1009 1005 13 ab fo <NA>
## 10 amp1010 1006 1 ab fo <NA>
## # … with 1,153 more rows, and 19 more variables: passive_methods <chr>,
## # complementary_methods <chr>, period <chr>, month_start <dbl>,
## # year_start <dbl>, month_finish <dbl>, year_finish <dbl>,
## # effort_months <dbl>, country <chr>, state <chr>, state_abbreviation <chr>,
## # municipality <chr>, site <chr>, latitude <dbl>, longitude <dbl>,
## # coordinate_precision <chr>, altitude <dbl>, temperature <dbl>,
## # precipitation <dbl>Por fim, podemos facilmente criar um objeto espacial do tipo MULTIPOINT utilizando a função sf::st_as_sf(). Podemos ver essas coordenadas plotadas no mapa simples da Figura @ref(fig:fig-vetor-pontos-atlantic-amphibians).
É necessário antes se ater primeiramente ao argumento coords que deve indicar as colunas de longitude e latitude, nessa ordem; e também ao argumento crs para indicar o CRS correspondente dessas coordendas, que aqui sabemos se tratar de coordenadas geográficas e datum WGS84. Então podemos facilmente utilizar o código EPSG 4326. Entretanto, se as coordenadas estiverem em metros, por exemplo, teremos de nos ater à qual CRS as mesmas foram coletadas, ou seja, se forem coordenadas de GPS, é preciso saber como o GPS estava configurado (projeção e datum).
# convert para sf
aa_locais_ve <- aa_locais %>%
sf::st_as_sf(coords = c("longitude", "latitude"), crs = 4326)
aa_locais_ve## Simple feature collection with 1163 features and 23 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -56.74194 ymin: -33.51083 xmax: -34.79667 ymax: -3.51525
## geographic CRS: WGS 84
## # A tibble: 1,163 x 24
## id reference_number species_number record sampled_habitat active_methods
## * <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 amp1001 1001 19 ab fo,ll as
## 2 amp1002 1002 16 co fo,la,ll as
## 3 amp1003 1002 14 co fo,la,ll as
## 4 amp1004 1002 13 co fo,la,ll as
## 5 amp1005 1003 30 co fo,ll,br as
## 6 amp1006 1004 42 co tp,pp,la,ll,is <NA>
## 7 amp1007 1005 23 co sp as
## 8 amp1008 1005 19 co sp,la,sw as,sb,tr
## 9 amp1009 1005 13 ab fo <NA>
## 10 amp1010 1006 1 ab fo <NA>
## # … with 1,153 more rows, and 18 more variables: passive_methods <chr>,
## # complementary_methods <chr>, period <chr>, month_start <dbl>,
## # year_start <dbl>, month_finish <dbl>, year_finish <dbl>,
## # effort_months <dbl>, country <chr>, state <chr>, state_abbreviation <chr>,
## # municipality <chr>, site <chr>, coordinate_precision <chr>, altitude <dbl>,
## # temperature <dbl>, precipitation <dbl>, geometry <POINT [°]>plot(aa_locais_ve[1], pch = 20, col = "black", main = NA, axes = TRUE, graticule = TRUE)
Coordenadas das comunidades do Atlantic Amphinians (Vancine et al. 2018).
Converter dados espaciais sp para sf
O pacote sf é mais recente e mais fácil de manipular objetos vetoriais no R, como vimos. Seu predecessor, o pacote sp possui uma classe própria e homônima. Entretanto, muitos pacotes de análises espaciais ainda utilizam essa classe em suas funções, apesar dessa migração ter ocorrido rapidamente recentemente. Dessa forma, a conversão entre essas classes pode ser necessária em alguns momentos.
Abaixo, veremos como podemos fazer essa conversão facilmente. Primeiramente, vamos importar dados sp.
# paises sp
co110_sp <- rnaturalearth::countries110
class(co110_sp)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"Agora, podemos converter facilmente com a função sf::st_as_sf().
# paises sf
co110_sf <- sf::st_as_sf(co110_sp)
class(co110_sf)## [1] "sf" "data.frame"Podemo facilmente converter esse objeto novamente para a classe sp com a função sf::as_Spatial.
# paises sp
co110_sp <- sf::as_Spatial(co110_sf)## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on WGS84 ellipsoid in Proj4
## definitionclass(co110_sp)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"Raster
Para importar dados raster no R, iremos utilizar a função raster::raster(), raster::brick() ou raster::stack(). Para apenas uma camada raster, usaremos a função raster::raster(), com o argumento x sendo o nome do arquivo. Já para mais camadas, usaremos raster::brick() para um arquivo que possua múltiplas camadas, ou ainda a função raster::stack() para várias arquivos em diferentes camadas também no argumento x, sendo necessário listar os arquivos no diretório, geralmente utilizando a função dir() ou list.files(). Entretanto, para especificar uma camada, podemos utiliar o argumento band ou layer e o nome dessa camada.
Raster Layer
Primeiramente, vamos criar um diretório como para os dados raster que iremos fazer o download.
# criar directorio
dir.create(here::here("dados", "raster"))Em seguida, vamos fazer o download de dados de elevação, na verdade dados de Modelo Digital de Elevação (Digital Elevation Model - DEM), localizados também para o município de Rio Claro. Iremos utilizar os dados do Shuttle Radar Topography Mission - SRTM. Para saber mais sobre esses dados, recomendamos a leitura do artigo Farr et al. (2007).
# aumentar o tempo de download
options(timeout = 600)
# download
download.file(url = "https://srtm.csi.cgiar.org/wp-content/uploads/files/srtm_5x5/TIFF/srtm_27_17.zip",
destfile = here::here("dados", "raster", "srtm_27_17.zip"), mode = "wb")
# unzip
unzip(zipfile = here::here("dados", "raster", "srtm_27_17.zip"),
exdir = here::here("dados", "raster"))Agora podemos importar essa camada para o R, e visualizá-la em relação ao limite do município de Rio Claro/SP (Figura @ref(fig:fig-raster-dem)).
# importar raster
ra <- raster::raster(here::here("dados", "raster", "srtm_27_17.tif"))
ra## class : RasterLayer
## dimensions : 6000, 6000, 3.6e+07 (nrow, ncol, ncell)
## resolution : 0.0008333333, 0.0008333333 (x, y)
## extent : -50, -45, -25, -20 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : /home/mude/Downloads/cap_14/dados/raster/srtm_27_17.tif
## names : srtm_27_17
## values : -32768, 32767 (min, max)plot(ra, col = viridis::viridis(10))
plot(rc_2019[1], col = NA, add = TRUE)
Camada raster do DEM em relação ao limite do município de Rio Claro/SP.
Raster Stack
Além dos dados de elevação, dados de temperatura e precipitação podem ser obtidos do WorldClim. Para saber mais sobre esses dados, recomandamos a leitura do artigo Fick & Hijmans (2017).
# download
download.file(url = "https://biogeo.ucdavis.edu/data/worldclim/v2.1/base/wc2.1_10m_bio.zip",
destfile = here::here("dados", "raster", "wc2.0_10m_bio.zip"), mode = "wb")
# unzip
unzip(zipfile = here::here("dados", "raster", "wc2.0_10m_bio.zip"),
exdir = here::here("dados", "raster"))Para importar essa série de camadas, primeiramente iremos listar os arquivos e depois importar no formato RasterStack ((Figura @ref(fig:fig-raster-wc)).
# listar arquivos
fi <- dir(path = here::here("dados", "raster"), pattern = "wc") %>%
grep(".tif", ., value = TRUE)
fi## [1] "wc2.1_10m_bio_1.tif" "wc2.1_10m_bio_10.tif" "wc2.1_10m_bio_11.tif"
## [4] "wc2.1_10m_bio_12.tif" "wc2.1_10m_bio_13.tif" "wc2.1_10m_bio_14.tif"
## [7] "wc2.1_10m_bio_15.tif" "wc2.1_10m_bio_16.tif" "wc2.1_10m_bio_17.tif"
## [10] "wc2.1_10m_bio_18.tif" "wc2.1_10m_bio_19.tif" "wc2.1_10m_bio_2.tif"
## [13] "wc2.1_10m_bio_3.tif" "wc2.1_10m_bio_4.tif" "wc2.1_10m_bio_5.tif"
## [16] "wc2.1_10m_bio_6.tif" "wc2.1_10m_bio_7.tif" "wc2.1_10m_bio_8.tif"
## [19] "wc2.1_10m_bio_9.tif"# importar
st <- raster::stack(here::here("dados", "raster", fi))
st## class : RasterStack
## dimensions : 1080, 2160, 2332800, 19 (nrow, ncol, ncell, nlayers)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : wc2.1_10m_bio_1, wc2.1_10m_bio_10, wc2.1_10m_bio_11, wc2.1_10m_bio_12, wc2.1_10m_bio_13, wc2.1_10m_bio_14, wc2.1_10m_bio_15, wc2.1_10m_bio_16, wc2.1_10m_bio_17, wc2.1_10m_bio_18, wc2.1_10m_bio_19, wc2.1_10m_bio_2, wc2.1_10m_bio_3, wc2.1_10m_bio_4, wc2.1_10m_bio_5, ...
## min values : -54.724354, -37.781418, -66.311249, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1.000000, 9.131122, 0.000000, -29.686001, ...
## max values : 30.98764, 38.21617, 29.15299, 11191.00000, 2381.00000, 484.00000, 229.00169, 5284.00000, 1507.00000, 5282.00000, 4467.00000, 21.14754, 100.00000, 2363.84595, 48.08275, ...plot(st[[c(1, 4)]], col = viridis::viridis(10))
Camadas rasters do WorldClim (BIO01 e BIO12) para o mundo.
Exportar dados
Saber a melhor forma de exportar dados geográficos de objetos recém-criados no R é fundamental, principalmente porque essa ação irá depender do tipo de dado (vetor ou raster), classe do objeto (por exemplo, MULTIPOINT ou RasterLayer) e tipo e quantidade de informações armazenadas (por exemplo, tamanho do objeto, intervalo de valores, etc.).
Vetor
Para dados vetoriais, a principal função utilizada é a sf::st_write(). Essa função permite gravar objetos sf em vários formatos de arquivos vetoriais, como .shp, .gpkg ou .geojson. O formato a ser exportado vai influenciar na velocidade do processo de gravação.
Os argumentos dessa função será o obj que é o objeto sf criado no ambiente R, e o dsn (data source name), ou seja, o nome que o arquivo terá ao ser exportado do R, de modo que o complemento .shp no nome de saída, por exemplo, irá definir que o arquivo terá a extensão ESRI Shapefile. Entretanto, essa extensão pode ser definida também utilizando o argumento driver, com as possibilidades listadas nesse site.
# exportar o vetor do limite de rio claro na extensão esri shapefile
sf::st_write(obj = rc_2019, dsn = here::here("dados", "vetor", "rio_claro.shp"))Ou podemos ainda exportar o objeto vetorial na extensão GeoPackage. Entretando, aqui é interessante acrescentar um argumento chamado layer para definir o nome das camadas a serem exportadas no mesmo arquivo GeoPackage, por exemplo.
# exportar o vetor do limite de rio claro na extensão geopackage
sf::st_write(obj = rc_2019, dsn = here::here("dados", "vetor", "vetores.gpkg"), layer = "rio_claro")Ainda sobre o formato GeoPackage, há algo muito interessante que podemos fazer: podemos agrescentar outros arquivos vetoriais ao mesmo arquivo já criado. Como exemplo, iremos exportar o limite do Brasil para o mesmo arquivo.
# exportar o vetor do limite do brasil na extensão esri shapefile
sf::st_write(obj = br, dsn = here::here("dados", "vetor", "vetores.gpkg"), layer = "brasil")Raster
Para exportar dados raster utilizamos geralmente a função raster::writeRaster(). Exportar dados raster é um pouco mais complexo que exportar dados vetoriais. Teremos de definir se iremos exportar arquivos em uma ou várias camadas, quantidade de informações por pixel, e ainda diferentes extensões de saída. Um ponto fundamental: arquivos raster escritos em discos geralmente ocupam bastante espaço, e dessa forma, há parâmetros específicos para certos tipos de dados, que detalharemos a seguir para contornar esse problema e comprimir os arquivos.
Na função raster::writeRaster(), o argumento x diz respeito ao objeto raster no ambiente R. O argumento filename é nome do arquivo que será exportado do R, podendo ou não possuir a extensão que se pretende que o arquivo tenha. O argumento format é o formato do arquivo, sendo as principais possibilidades resumidas na Tabela @ref(tab:tab-raster-formatos), e para saber das possibilidade suportadas, use a função raster::writeFormats(). O argumento bylayer diz se de um objeto com múltiplas camadas, cada uma delas será exportada em um arquivo diferente.
| Tipo de arquivo | Nome longo | Extensão | Suporte a múltiplas camadas |
|---|---|---|---|
| raster | Formato pacote raster | .grd | Sim |
| ascii | ESRI Ascii | .asc | Não |
| SAGA | SAGA GIS | .sdat | Não |
| IDRISI | IDRISI | .rst | Não |
| CDF | netCDF (requer ncdf4) | .nc | Sim |
| GTiff | GeoTiff (requer rgdal) | .tif | Sim |
| ENVI | ENVI .hdr | .envi | Sim |
| EHdr | ESRI .hdr | .bil | Sim |
| HFA | Erdas imagem (.img) | .img | Sim |
Dentre os argumentos adicionais, temos ainda o datatype, que faz referência à um dos nove tipos de dados detalhados na Tabela @ref(tab:tab-raster-tipos), sendo que o tipo de dado determina a representação em bits (quantidade de informação) na célula do objeto raster exportado e depende da faixa de valores do objeto raster em cada pixel. Quanto mais valores um tipo de dado puder representar, maior será o arquivo exportado no disco, dessa forma, é interessante utilizar um tipo de dado que diminua o tamanho do arquivo à ser exportado, dependendo do tipo de dado em cada pixel. Para a função raster::writeRaster(), o default é FLT4S, o que pode ocupar mais espaço em disco do que o necessário.
| Tipo de dado | Valor mínimo | Valor máximo |
|---|---|---|
| LOG1S | FALSE (0) | TRUE (1) |
| INT1S | -127 | 127 |
| INT1U | 0 | 255 |
| INT2S | -32.767 | 32.767 |
| INT2U | 0 | 65534 |
| INT4S | -2.147.483.647 | 2.147.483.647 |
| INT4U | 0 | 42.94.967.296 |
| FLT4S | -3,4e+38 | 3,4e+38 |
| FLT8S | -1,7e+308 | 1,7e+308 |
Outros argumentos de suporte são: overwrite para sobreescrever um arquivo que já exista, progress para mostrar uma barra de progresso da exportação como “text” ou “window”, e options que permite opções do GDAL. Para esse último argumento, quando exportar especificamente na extensão GeoTIFF, podemos utilizar options = c("COMPRESS=NONE", "TFW=YES") para que haja compressão do arquivo, diminuindo consideravelmente seu tamanho (cerca de um terço), aliado à criação de um arquivo auxiliar .tfw, para ser carregado em softwares específicos de SIG, como o ArcGIS.
Para exportar apenas uma camada RasterLayer, podemos utilizar a função raster::writeRaster() em um formato mais simples.
# diretorio
dir.create(here::here("dados", "raster", "exportados"))
# exportar raster layer
raster::writeRaster(ra,
filename = here::here("dados", "raster", "exportados", "elevation"),
format = "GTiff",
datatype = "INT2S",
options = c("COMPRESS=NONE", "TFW=YES"),
progress = "text",
overwrite = TRUE)Para mais de uma camada RasterBrick ou RasterStack, podemos utilizar a função raster::writeRaster() com mais argumentos, como o bylayer = TRUE.
raster::writeRaster(x = st,
filename = here::here("dados", "raster", "exportados", names(st)),
bylayer = TRUE,
format = "GTiff",
datatype = "INT2S",
options = c("COMPRESS=NONE", "TFW=YES"),
progress = "text",
overwrite = TRUE)Descrição de objetos espaciais
Muitas vezes iremos precisar verificar as informações dos objetos geográficos importados para o R. Apesar de chamar o objeto trazer grande parte das informações que precisamos consultar, existem funções específicas que nos auxiliam nesse processo de descrição dos objetos.
Vetor
Podemos acessar as informações geográficas e a tabela de atributos de um objeto importado como vetor simplesmente chamando o nome do objeto no R.
# rio claro
rc_2019## Simple feature collection with 1 feature and 7 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -47.76536 ymin: -22.55203 xmax: -47.46188 ymax: -22.24368
## geographic CRS: SIRGAS 2000
## code_muni name_muni code_state abbrev_state name_state code_region
## 493 3543907 Rio Claro 35 SP São Paulo 3
## name_region geom
## 493 Sudeste MULTIPOLYGON (((-47.66303 -...Mas também podemos acessar informações geográficas com funções específicas, como tipo de geometria, limites geográficos do vetor (extensão), sistema de referência de coordenadas (CRS), e a tabela de atributos.
# tipo de geometria
sf::st_geometry_type(rc_2019)## [1] MULTIPOLYGON
## 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE# extensao
sf::st_bbox(rc_2019)## xmin ymin xmax ymax
## -47.76536 -22.55203 -47.46188 -22.24368# crs
sf::st_crs(rc_2019)## Coordinate Reference System:
## User input: SIRGAS 2000
## wkt:
## GEOGCRS["SIRGAS 2000",
## DATUM["Sistema de Referencia Geocentrico para las AmericaS 2000",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["Latin America - Central America and South America - onshore and offshore. Brazil - onshore and offshore."],
## BBOX[-59.87,-122.19,32.72,-25.28]],
## ID["EPSG",4674]]# acessar a tabela de atributos
rc_2019_tab <- sf::st_drop_geometry(rc_2019)
rc_2019_tab## code_muni name_muni code_state abbrev_state name_state code_region
## 493 3543907 Rio Claro 35 SP São Paulo 3
## name_region
## 493 SudesteRaster
Da mesma forma, podemos acessar as informações objetos raster chamando o nome do objeto.
ra## class : RasterLayer
## dimensions : 6000, 6000, 3.6e+07 (nrow, ncol, ncell)
## resolution : 0.0008333333, 0.0008333333 (x, y)
## extent : -50, -45, -25, -20 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : /home/mude/Downloads/cap_14/dados/raster/srtm_27_17.tif
## names : srtm_27_17
## values : -32768, 32767 (min, max)Além disso, podemos selecionar informações desse objeto com funções específicas, tanto para RasterLayer, quanto para RasterBrick ou RasterStack como: classe, dimensões (número de linhas, colunas e camadas), número de camadas, número de linhas, número de colunas, número de células, resolução (largura e altura do tamanho do pixel), extensão (limites geográficos), sistema de referência de coordenadas (CRS), nome das camadas e extrair os valores de todos os pixels.
# classe
class(ra)## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"# dimensoes
dim(ra)## [1] 6000 6000 1# numero de camadas
nlayers(ra)## [1] 1# numero de linhas
nrow(ra)## [1] 6000# numero de colunas
ncol(ra)## [1] 6000# numero de celulas
ncell(ra)## [1] 3.6e+07# resolucao
res(ra)## [1] 0.0008333333 0.0008333333# extensao
extent(ra)## class : Extent
## xmin : -50
## xmax : -45
## ymin : -25
## ymax : -20# projecao ou crs
projection(ra)## [1] "+proj=longlat +datum=WGS84 +no_defs"# nome
names(ra)## [1] "srtm_27_17"# valores
getValues(ra) %>% head## [1] 382 379 379 379 379 383values(ra) %>% head## [1] 382 379 379 379 379 383ra[] %>% head## [1] 382 379 379 379 379 383Reprojeção de dados geográficos
Em algumas situações é necessário alterar o CRS de um objeto espacial para um novo CRS. A reprojeção é justamente a transformação de coordenadas de um CRS para outro: geográficos (‘lon/lat’, com unidades em graus de longitude e latitude) e projetados (normalmente com unidades de metros a partir de um datum).
Geralmente iremos precisar fazer essa operação para transformar camadas vetoriais ou rasters para o mesmo CRS, de modo que possam ser exibidas conjuntamente, ou ainda que as camadas possuem CRS projetado para realizar alguma operação espacial entre camadas, ou quando precisamos calcular áreas, formatos ou distâncias, como métricas de paisagem, por exemplo. Existe uma infinidade de projeções e um excelente material de consulta é o livro de Lapaine et al. (2017).
Podemos verificar o CRS de uma camada através da função sf::st_crs() ou raster::projection() e raster::crs(), ou ainda, saber se a mesma possui um CRS geográfico ou não, com a função sf::st_is_longlat().
Já para reprojetar um objeto sf usamos a função sf::st_transform() e para um objeto raster usamos a função raster::projectRaster().
# projecao de vetores
sf::st_crs(rc_2019)## Coordinate Reference System:
## User input: SIRGAS 2000
## wkt:
## GEOGCRS["SIRGAS 2000",
## DATUM["Sistema de Referencia Geocentrico para las AmericaS 2000",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["Latin America - Central America and South America - onshore and offshore. Brazil - onshore and offshore."],
## BBOX[-59.87,-122.19,32.72,-25.28]],
## ID["EPSG",4674]]# projecao de raster
raster::projection(ra)## [1] "+proj=longlat +datum=WGS84 +no_defs"raster::crs(ra)## CRS arguments: +proj=longlat +datum=WGS84 +no_defs# verificar se o crs e geografico
sf::st_is_longlat(rc_2019)## [1] TRUEAs funções sf::st_transform() e raster::projectRaster() possuem dois parâmetros importantes: x que é o objeto a ser reprojetado e o crs que é o CRS alvo. O argumento crs pode ser especificado de quatro maneiras: 1) código EPSG (por exemplo, 4326), 2) string PROJ4 (por exemplo, “+ proj = longlat + datum = WGS84 + no_defs”), 3) string WKT, ou 4) objeto crs de outra camada, conforme retornado por sf::st_crs() ou raster::crs(). Esas informações de EPSG, PROJ4 e WKT podem ser acessadas nas bases: epsg.io e spatialreference.org.
Dentre os possíveis CRSs a serem utilizados, alguns são mais comuns para CRSs geográficos e projetados. Para CRSs geográficos, o mais comum para o mundo é o World Geodetic System 1984 (WGS84), ou seja, geográfico com datum WGS84. Para o Brasil, o CRS adotado é o Sistema de Referencia Geocéntrico para las Américas 2000 (SIRGAS 2000), ou seja, geográfico com datum SIRGAS2000.
Para CRSs projetados, essa escolha vai depender da extenção e localização da área de interesse no globo terrestre. Aqui destacaremos os principais, para três escalas: global, regional e local. Para a escala global, geralmente usa-se umas dessas projeções, dependendo do objetivo: 1) Projeção de Mollweide, 2) Projeção de Winkel Tripel, 3) Projeção de Eckert IV, 4) Projeção Azimutal de Lambert. Para a escala regional, como um hemisfério, geralmente usa-se a Projeção Cônica de Albers. Por fim, para a escala local, usa-se geralmente a Projeção Universal Transverse Mercator (UTM), um conjunto de CRSs que divide a Terra em 60 cunhas longitudinais e 20 segmentos latitudinais, como pode ser visto neste link.
Os principais CRSs são descritos na Tabela @ref(tab:tab-crs).
| CRS | Tipo de CRS | Descrição | epsg.io | spatialreference.org |
|---|---|---|---|---|
| World Geodetic System 1984 (WGS84) | Geográfico | CRS geográfico mais comum para o mundo | EPSG:4326 | EPSG:4326 |
| Sistema de Referencia Geocéntrico para las Américas 2000 (SIRGAS 2000) | Geográfico | CRS geográfico oficial para o Brasil | EPSG:4674 | EPSG:4674 |
| Projeção de Mollweide | Projetado | CRS projetado que preserva as relações de área | ESRI:54009 | SR-ORG:7099 |
| Projeção de Winkel Tripel | Projetado | CRS projetado com mínimo de distorção para área, direção e distância | NA | SR-ORG:7291 |
| Projeção de Eckert IV | Projetado | CRS projetado que presenva a área e com meridianos elípticos | EPSG:54012 | ESRI:54012 |
| Projeção Azimutal de Lambert | Projetado | CRS projetado que preserva os tamanhos relativos e senso de direção a partir do centro | NA | NA |
| Projeção Cônica de Albers | Projetado | CRS projetado para escala regional, mantendo a área constante em toda sua superfície | NA | SR-ORG:7823 |
| Projeção Universal Transverse Mercator (UTM) | Projetado | CRS projetado para escala local, distorcendo áreas e distâncias com gravidade crescente com a distância do centro da zona UTM | EPSG:31983 | EPSG:31983 |
Vetor
Como dissemos anteriormente, para reprojetar um vetor, utilizamos a função sf::st_transform(), observando os argumentos x que é a camada a ser reprojetada, e o crs que é o CRS alvo.
Vamos reprojetar o limite do município de Rio Claro/SP do CRS SIRGAS2000/geográfico para o CRS projetado SIRGAS2000/UTM23S, com os efeitos da transformação podendo ser notados na Figura @ref(fig:fig-vetor-crs-trans).
# converter crs
rc_2019_sirgas2000_utm23s <- sf::st_transform(x = rc_2019, crs = 31983)
Limites do município de Rio Claro/SP com CRS SIRGAS2000/geográfico e com CRS SIRGAS2000/UTM23S.
Podemos ainda utilizar o formato proj4string no argumento crs para fazer a transformação. Vamos primeiramente plotar o mundo em WGS84/Geográfico (Figura @ref(fig:fig-vetor-mundo-wgs84)).
plot(co110_sf[1], col = "gray", main = "WGS84/Geográfio", graticule = TRUE)
Camada BIO01 para o mundo com CRS geográfico e datum WGS84.
Agora, iremos reprojetar utilizando a Projeção de Mollweide (Figura @ref(fig:fig-vetor-mundo-moll)).
# projecao de mollweide
co110_sf_moll <- sf::st_transform(x = co110_sf, crs = "+proj=moll")plot(co110_sf_moll[1], col = "gray", main = "Projeção de Mollweide", graticule = TRUE)
Camada BIO01 para o mundo com CRS Projeção de Mollweide.
Ou ainda podemos utilizar a Projeção Azimutal de Lambert com alguns parâmetros ajustados para centralizar a projeção no Brasil (@ref(fig:fig-vetor-mundo-laea)).
# projecao de mollweide
co110_sf_laea <- sf::st_transform(x = co110_sf,
crs = "+proj=laea +x_0=0 +y_0=0 +lon_0=-50 +lat_0=0")plot(co110_sf_laea[1], col = "gray", main = "Projeção Azimutal de Lambert", graticule = TRUE)
amada BIO01 para o mundo com CRS Projeção Azimutal de Lambert centrado no Brasil.
Raster
A reprojeção de objetos raster não é uma tarefa tão simples quanto a reprojeção de vetores. Em vetores, a reprojeção altera as coordenadas de cada vértice. Entretanto, como rasters são compostos de células retangulares do mesmo tamanho, a reprojeção do raster envolve a criação de um novo objeto raster, envolvendo duas operações espaciais separadas: 1) reprojeção vetorial dos centróides celulares para outro CRS (i.e., muda a posição e tamanho do pixel) e, 2) cálculo de novos valores do pixel por meio de reamostragem (i.e., muda o valor do pixel).
A função raster::projectRaster() possui alguns parâmetros que necessitam de algumas especificações. O argumento from que é o objeto raster de entrada que sofre a reprojeção. O argumento to é um objeto raster do qual todas as propriedade CRSs, como extenção e resolução serão associadas ao objeto raster indicado em from. O argumento res permite ajustar a resolução do pixel de saída do objeto raster reprojetado.
O argumento crs aceita apenas as definições de proj4string extensas de um CRS em vez de códigos EPSG concisos. Contudo, é possível usar um código EPSG em uma definição de proj4string com +init=epsg:EPSG. Por exemplo, pode-se usar a definição +init=epsg:4326 para definir CRS para WGS84 (código EPSG de 4326). A biblioteca PROJ adiciona automaticamente o resto dos parâmetros e os converte em +init=epsg:4326 +proj=longlat +datum=WGS84 + no_defs + ellps=WGS84 + towgs84=0,0,0.
O argumento method permite escolher entre os métodos “ngb” (vizinho mais próximo) ou “biliniar” (interpolação bilinear), sendo o primeiro mais indicado para reprojeção de rasters categóricos, pois os valores estimados devem ser iguais aos do raster original. O método “ngb” define cada novo valor de célula para o valor da célula mais próxima (centro) do raster de entrada. Já o método “biliniar” é indicado para raster contínuos e calcula o valor da célula de saída com base nas quatro células mais próximas no raster original, sendo a média ponderada da distância dos valores dessas quatro células.
Aqui, vamos reprojetar os dados de elevação para Rio Claro/SP. Para que esse processo seja mais rápido, iremos antes ajustar a extensão do raster para o limite do município usando a função raster::crop() (Figura @ref(fig:fig-raster-crop)). Essa função é melhor explicada na seção xx.
# ajuste do limite
ra_rc <- raster::crop(x = ra, y = rc_2019)
ra_rc## class : RasterLayer
## dimensions : 370, 364, 134680 (nrow, ncol, ncell)
## resolution : 0.0008333333, 0.0008333333 (x, y)
## extent : -47.765, -47.46167, -22.55167, -22.24333 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : srtm_27_17
## values : 491, 985 (min, max)plot(ra_rc, col = viridis::viridis(10))
plot(rc_2019[1], col = NA, lwd = 2, add = TRUE)
Ajuste da extensão do raster de elevação para o município de Rio Claro/SP.
Primeiramente, vamos reprojetar indicando uma projeção e sem especificar o tamanho da célula. Note que o tamanho da célula vai se ajustar para valores diferentes, sendo portanto, pixels retangulares e não quadrados.
# reprojecao
ra_rc_sirgas2000_utm23s <- raster::projectRaster(from = ra_rc, crs = "+init=epsg:31983", method = "bilinear")## Warning in showSRID(uprojargs, format = "PROJ", multiline
## = "NO", prefer_proj = prefer_proj): Discarded datum
## Sistema_de_Referencia_Geocentrico_para_las_AmericaS_2000 in Proj4 definitionra_rc_sirgas2000_utm23s## class : RasterLayer
## dimensions : 386, 381, 147066 (nrow, ncol, ncell)
## resolution : 85.8, 92.3 (x, y)
## extent : 214575.4, 247265.2, 7503009, 7538637 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=23 +south +ellps=GRS80 +units=m +no_defs
## source : memory
## names : srtm_27_17
## values : 491.6033, 980.4151 (min, max)Agora vamos reprojetar especificando o tamanho da célula (Figura @ref(fig:fig-raster-reproj)). Dessa forma, todas as células terão o mesmo, i.e., quadrados de 90 metros.
# reprojecao
ra_rc_sirgas2000_utm23s <- raster::projectRaster(from = ra_rc, crs = "+init=epsg:31983", method = "bilinear", res = 90)## Warning in showSRID(uprojargs, format = "PROJ", multiline
## = "NO", prefer_proj = prefer_proj): Discarded datum
## Sistema_de_Referencia_Geocentrico_para_las_AmericaS_2000 in Proj4 definitionra_rc_sirgas2000_utm23s## class : RasterLayer
## dimensions : 396, 364, 144144 (nrow, ncol, ncell)
## resolution : 90, 90 (x, y)
## extent : 214554.4, 247314.4, 7502985, 7538625 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=23 +south +ellps=GRS80 +units=m +no_defs
## source : memory
## names : srtm_27_17
## values : 493.2395, 986.686 (min, max)plot(ra_rc_sirgas2000_utm23s, col = viridis::viridis(10))
plot(rc_2019_sirgas2000_utm23s[1], col = NA, lwd = 2, add = TRUE)
Reprojeção do raster de elevação para SIRGAS2000/UTM23S especificado por um objeto e informando o tamanho da célula.
Vamos também reprojetar uma camada mundial da média de temperatura anual (BIO01), indicando o tamanho da célula para 25.000 m (Figura @ref(fig:fig-raster-reproj-celula-mundo)).
# reprojecao
bio01_moll <- raster::projectRaster(st[[1]], crs = "+proj=moll", res = 25000, method = "bilinear")## Warning in rgdal::rawTransform(projto_int, projfrom, nrow(xy), xy[, 1], : 245665
## projected point(s) not finitebio01_moll## class : RasterLayer
## dimensions : 732, 1453, 1063596 (nrow, ncol, ncell)
## resolution : 25000, 25000 (x, y)
## extent : -18159905, 18165095, -9154952, 9145048 (xmin, xmax, ymin, ymax)
## crs : +proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## source : memory
## names : wc2.1_10m_bio_1
## values : -54.66752, 30.71805 (min, max)plot(bio01_moll, col = viridis::viridis(10))
plot(co110_sf_moll[1], col = NA, add = TRUE)
Reprojeção do raster de média de temperatura anual (BIO01) para Projeção de Mollweide informando o tamanho da célula.
Principais operações com dados geográficos
Nesta seção veremos as principais funções para realizar operações com dados geográficos. Essas operações são separadas conforme Lovelace, Nowosad & Muenchow (2019) em: Operações de atributos, Operações espaciais, e Operações geométricas.
Operações de atributos
São modificação de objetos espaciais baseado em informações não espaciais associadas a dados geográficos, como a tabela de atributos ou valores das células e nome das camadas dos rasters.
Vetor
As principais operações de atributos vetoriais são com respeito à tabela de atributos, sendo elas: 1) filtro, 2) junção, 3) agregação, e 4) manipulação da tabela de atributos. A lista de possíveis operações é longa, dessa forma, apresentaremos algumas operações utilizando as princiais funções e listamos as demais funções e suas operações, que irão depender de objetivos específicos.
Quase todas as operações serão as mesmas realizadas pelo pacote dplyr em uma tabela de dados (ver serção xx), sendo algumas operações específicas para alterar apenas campos da tabela de atributos e outras que refletem operações nas feições, ou seja, irão alterar através da tabela de atributos as características das feições. Essas funções e suas operações são descritas com detalhes na Tabela (@ref(tab:tab-vetor-operacoes-atributos)).
| Funções | Onde atua | Descrição |
|---|---|---|
filter() |
Feições | Selecionar feições por valores |
slice() |
Feições | Selecionar feições pela posição na tabela de atributos |
n_sample() |
Feições | Amostrar feições na tabela de atributos |
group_by() |
Feições | Agrupar feições por valores da tabela de atributos |
summarise() |
Feições | Operações com valores das feições na tabela de atributos, que acabam por dissolver as feições |
select() |
Atributos | Selecionar colunas da tabela de atributos |
pull() |
Atributos | Selecionar uma coluna da tabela de atributos como vetor |
rename() |
Atributos | Renomear uma coluna da tabela de atributos |
mutate() |
Atributos | Criar uma coluna ou alterar os valores da tabela de atributos |
*_join() |
Atributos | Diversas funções para juntar dados de outras tabelas de dados à tabela de atributos |
Para exemplificar as operações de atributos, vamos utilizar os dados de nascentes, hidrologia e cobertura da terra para o município de Rio Claro/SP.
Filtro
Vamos iniciar as operações fazendo o filtro de feições pela tabela de atributos, que permite selecionar feições pelos seus valores atribuídos, utilizando a função dplyr::filter(). Aqui vamos selecionar as feições de floresta do mapa de cobertura da terra para Rio Claro/SP (Figura @ref(fig:fig-vetor-opat-filtro)).
# filtro
rc_cob_floresta <- rc_cob %>%
dplyr::filter(CLASSE_USO == "formação florestal")plot(rc_2019_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(rc_cob_floresta$geometry, col = "forestgreen", add = TRUE)
Filtro da classe floresta para o mapeamento de cobertura da terra para o município de Rio Claro/SP.
Junção
Uma das funções mais úteis de operações de atributos é a junção, referida em inglês como join, através das funções dplyr::*_join() (ver detalhes na seção xx). Nela, usamos uma coluna identificadora para atribuir dados de outra tabela de dados. Como exemplo, vamos criar uma tabela de dados com novos nomes das classes de cobertura da terra e atribuir esses novos nomes à tabela de atributos do objeto vetorial. É fundamental destacar que para que essa função funcione, precisamos de uma coluna identificadora dos valores para que a junção seja possível.
# dados
da_classes <- tibble::tibble(CLASSE_USO = rc_cob$CLASSE_USO,
classe = c("agua", "antropico", "edificado", "floresta", "silvicultura"))
da_classes## # A tibble: 5 x 2
## CLASSE_USO classe
## <chr> <chr>
## 1 água agua
## 2 área antropizada antropico
## 3 área edificada edificado
## 4 formação florestal floresta
## 5 silvicultura silvicultura# juncao
rc_cob_classes <- dplyr::left_join(rc_cob, da_classes, by = "CLASSE_USO") %>%
sf::st_drop_geometry()
rc_cob_classes## GEOCODIGO MUNICIPIO UF CD_UF CLASSE_USO AREA_HA classe
## 1 3543907 RIO CLARO SP 35 água 357.027 agua
## 2 3543907 RIO CLARO SP 35 área antropizada 37297.800 antropico
## 3 3543907 RIO CLARO SP 35 área edificada 5078.330 edificado
## 4 3543907 RIO CLARO SP 35 formação florestal 7017.990 floresta
## 5 3543907 RIO CLARO SP 35 silvicultura 138.173 silviculturaAgregação
Outra função bastante útil é a agregação de atributos. Apesar de existir uma função que realiza a união de feições, o uso conjunto das funções dplyr::group_by() e dplyr::summarise() realizam uma tarefa semelhante. Aqui vamos agregar as nascentes para Rio Claro/SP, i.e., juntar cada ponto que estava numa linha da tabela de atributos de modo que todos fiquem numa mesma linha, com o valor da quantidade de nascentes (Figura @ref(fig:fig-vetor-opat-agregar)).
# agregar
rc_nas_n <- rc_nas %>%
dplyr::group_by(MUNICIPIO, HIDRO) %>%
dplyr::summarise(n = n())## `summarise()` has grouped output by 'MUNICIPIO'. You can override using the `.groups` argument.rc_nas_n## Simple feature collection with 1 feature and 3 fields
## geometry type: MULTIPOINT
## dimension: XY
## bbox: xmin: 217622.9 ymin: 7504132 xmax: 246367.4 ymax: 7537855
## projected CRS: SIRGAS 2000 / UTM zone 23S
## # A tibble: 1 x 4
## # Groups: MUNICIPIO [1]
## MUNICIPIO HIDRO n geometry
## <chr> <chr> <int> <MULTIPOINT [m]>
## 1 RIO CLARO nascen… 1220 ((217622.9 7528315), (217836.5 7528103), (217988.9 75…plot(rc_2019_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(rc_nas_n$geometry, pch = 20, col = "blue", add = TRUE)
Agregação e contagem das nascentes para o município de Rio Claro/SP.
Manipulação da tabela de atributos
Por fim, é muito comum em análises de softwares SIG a criação ou atualização dos valores de colunas na tabela de atributos. Aqui, podemos utilizar a função dplyr::mutate() para criar essas novas colunas, assim como atualizar os valores de colunas existentes. Em nosso exemplo, iremos fazer uma composição das colunas CLASSE_USO e AREA_HA numa nova coluna chamada classe_area.
# criar coluna
rc_cob_cob_col_area <- rc_cob %>%
dplyr::mutate(classe_area = paste0(CLASSE_USO, " (", AREA_HA, " ha)")) %>%
sf::st_drop_geometry()
rc_cob_cob_col_area## GEOCODIGO MUNICIPIO UF CD_UF CLASSE_USO AREA_HA
## 1 3543907 RIO CLARO SP 35 água 357.027
## 2 3543907 RIO CLARO SP 35 área antropizada 37297.800
## 3 3543907 RIO CLARO SP 35 área edificada 5078.330
## 4 3543907 RIO CLARO SP 35 formação florestal 7017.990
## 5 3543907 RIO CLARO SP 35 silvicultura 138.173
## classe_area
## 1 água (357.027 ha)
## 2 área antropizada (37297.8 ha)
## 3 área edificada (5078.33 ha)
## 4 formação florestal (7017.99 ha)
## 5 silvicultura (138.173 ha)Duas funções são bastante interessantes de serem integradas junto a manipulação de tabelas de atributos. Elas calculam propriedades geométricas numéricas dos vetores de linhas (comprimento) e polígonos (área): sf::st_length() e sf::st_area(). Essas funções calculam essas propriedades em metros para comprimento e em metros quadrados para área, independentemente do CRS. Para tanto, vamos utilizar as linhas de hidrografia e os polígonos de cobertura da terra para Rio Claro/SP, e atribuir esses valores à tabela de atributos de ambos os objetos espaciais, utilizando em conjunto a função dplyr::mutate().
# comprimento das linhas
rc_hid_comp <- rc_hid %>%
dplyr::mutate(comprimento = sf::st_length(.))
rc_hid_comp## Simple feature collection with 1 feature and 7 fields
## geometry type: MULTILINESTRING
## dimension: XY
## bbox: xmin: 215155.3 ymin: 7504132 xmax: 246367.4 ymax: 7537978
## projected CRS: SIRGAS 2000 / UTM zone 23S
## GEOCODIGO MUNICIPIO UF CD_UF HIDRO COMP_KM
## 1 3543907 RIO CLARO SP 35 curso d'água (0 - 10m) 1142.98
## geometry comprimento
## 1 MULTILINESTRING ((231815.7 ... 1142981 [m]# area das poligonos
rc_cob_area <- rc_cob %>%
dplyr::mutate(area_m2 = sf::st_area(.))
rc_cob_area## Simple feature collection with 5 features and 7 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 215151.7 ymin: 7503723 xmax: 246582.4 ymax: 7537978
## projected CRS: SIRGAS 2000 / UTM zone 23S
## GEOCODIGO MUNICIPIO UF CD_UF CLASSE_USO AREA_HA
## 1 3543907 RIO CLARO SP 35 água 357.027
## 2 3543907 RIO CLARO SP 35 área antropizada 37297.800
## 3 3543907 RIO CLARO SP 35 área edificada 5078.330
## 4 3543907 RIO CLARO SP 35 formação florestal 7017.990
## 5 3543907 RIO CLARO SP 35 silvicultura 138.173
## geometry area_m2
## 1 MULTIPOLYGON (((235487.6 75... 3570267 [m^2]
## 2 MULTIPOLYGON (((232275 7504... 372978415 [m^2]
## 3 MULTIPOLYGON (((233123.6 75... 50783283 [m^2]
## 4 MULTIPOLYGON (((232355 7504... 70179895 [m^2]
## 5 MULTIPOLYGON (((243052.1 75... 1381726 [m^2]Raster
Devido a estrutura espacial do rater ser formada por uma ou mais superfícies contínuas, as manipulações como subconjunto e outras operações em objetos raster funcionam de uma maneira diferente do que em objetos vetoriais. Veremos aqui as três principais: 1) subconjunto de células usando o operador [] ou subconjunto de camadas RasterStack ou RasterBrick utilizando os operadores [[]] e $, 2) renomear nomes das camdas, e 3) resumir informações de todos os pixels.
Subconjunto
Podemos fazer um subconjunto de células utilizando dentro dos operadores [] valores para indicar a posição da linha e da coluna de um raster, ou ainda a posição de uma célula utilizando apenas um número. Essas operações resultarão em valores diferentes para RasterLayer e RasterBrick ou RasterStack.
# raster - linha 1 e columna 1
ra[1, 1]##
## 382# celula 1
ra[1]##
## 382# stack - linha 1 e columna 1
st[1, 1]## wc2.1_10m_bio_1 wc2.1_10m_bio_10 wc2.1_10m_bio_11 wc2.1_10m_bio_12
## [1,] NA NA NA NA
## wc2.1_10m_bio_13 wc2.1_10m_bio_14 wc2.1_10m_bio_15 wc2.1_10m_bio_16
## [1,] NA NA NA NA
## wc2.1_10m_bio_17 wc2.1_10m_bio_18 wc2.1_10m_bio_19 wc2.1_10m_bio_2
## [1,] NA NA NA NA
## wc2.1_10m_bio_3 wc2.1_10m_bio_4 wc2.1_10m_bio_5 wc2.1_10m_bio_6
## [1,] NA NA NA NA
## wc2.1_10m_bio_7 wc2.1_10m_bio_8 wc2.1_10m_bio_9
## [1,] NA NA NA# celula 1
st[1]## wc2.1_10m_bio_1 wc2.1_10m_bio_10 wc2.1_10m_bio_11 wc2.1_10m_bio_12
## [1,] NA NA NA NA
## wc2.1_10m_bio_13 wc2.1_10m_bio_14 wc2.1_10m_bio_15 wc2.1_10m_bio_16
## [1,] NA NA NA NA
## wc2.1_10m_bio_17 wc2.1_10m_bio_18 wc2.1_10m_bio_19 wc2.1_10m_bio_2
## [1,] NA NA NA NA
## wc2.1_10m_bio_3 wc2.1_10m_bio_4 wc2.1_10m_bio_5 wc2.1_10m_bio_6
## [1,] NA NA NA NA
## wc2.1_10m_bio_7 wc2.1_10m_bio_8 wc2.1_10m_bio_9
## [1,] NA NA NAPara selecionar uma camada de um RasterBrick ou RasterStack, podemos utilizar as funções raster::subset() ou raster::raster() com o argumento layer indicando a ordem ou o nome da camada, além dos operadores [[]] e $ (Figura @ref(fig:fig-raster-stack-subset)).
# selecao de camada num objeto stack utilizando a funcao subset
st_bio01 <- raster::subset(st, "wc2.1_10m_bio_1")
st_bio01## class : RasterLayer
## dimensions : 1080, 2160, 2332800 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : /home/mude/Downloads/cap_14/dados/raster/wc2.1_10m_bio_1.tif
## names : wc2.1_10m_bio_1
## values : -54.72435, 30.98764 (min, max)# selecao de camada num objeto stack utilizando a funcao raster
st_bio01 <- raster::raster(st, layer = 1)
st_bio01## class : RasterLayer
## dimensions : 1080, 2160, 2332800 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : /home/mude/Downloads/cap_14/dados/raster/wc2.1_10m_bio_1.tif
## names : wc2.1_10m_bio_1
## values : -54.72435, 30.98764 (min, max)# selecao de camada num objeto stack utilizando os operadores [[]] e o nome
st_bio01 <- st[["wc2.1_10m_bio_1"]]
st_bio01## class : RasterLayer
## dimensions : 1080, 2160, 2332800 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : /home/mude/Downloads/cap_14/dados/raster/wc2.1_10m_bio_1.tif
## names : wc2.1_10m_bio_1
## values : -54.72435, 30.98764 (min, max)# selecao de camada num objeto stack utilizando os operadores [[]] e a posicao
st_bio01 <- st[[1]]
st_bio01## class : RasterLayer
## dimensions : 1080, 2160, 2332800 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : /home/mude/Downloads/cap_14/dados/raster/wc2.1_10m_bio_1.tif
## names : wc2.1_10m_bio_1
## values : -54.72435, 30.98764 (min, max)# selecao de camada num objeto stack utilizando o operador $
st_bio01 <- st$wc2.1_10m_bio_1
st_bio01## class : RasterLayer
## dimensions : 1080, 2160, 2332800 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : /home/mude/Downloads/cap_14/dados/raster/wc2.1_10m_bio_1.tif
## names : wc2.1_10m_bio_1
## values : -54.72435, 30.98764 (min, max)raster::plot(st_bio01, col = viridis::viridis(10))
Camada BIO01 selecionada pelas operações de subconjunto acima do stack de variáveis bioclimáticas.
Renomear
Podemos ainda renomear camadas dos raster RasterLayer utilizando a função names().
# nomes
names(ra_rc)## [1] "srtm_27_17"# renomear
names(ra_rc) <- "elevacao"
# nomes
names(ra_rc)## [1] "elevacao"E essa operação também funciona para RasterBrick e RasterStack.
# nomes
names(st)## [1] "wc2.1_10m_bio_1" "wc2.1_10m_bio_10" "wc2.1_10m_bio_11" "wc2.1_10m_bio_12"
## [5] "wc2.1_10m_bio_13" "wc2.1_10m_bio_14" "wc2.1_10m_bio_15" "wc2.1_10m_bio_16"
## [9] "wc2.1_10m_bio_17" "wc2.1_10m_bio_18" "wc2.1_10m_bio_19" "wc2.1_10m_bio_2"
## [13] "wc2.1_10m_bio_3" "wc2.1_10m_bio_4" "wc2.1_10m_bio_5" "wc2.1_10m_bio_6"
## [17] "wc2.1_10m_bio_7" "wc2.1_10m_bio_8" "wc2.1_10m_bio_9"# renomear
names(st) <- c("bio01", paste0("bio", 10:19), paste0("bio0", 2:9))
# nomes
names(st)## [1] "bio01" "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16" "bio17"
## [10] "bio18" "bio19" "bio02" "bio03" "bio04" "bio05" "bio06" "bio07" "bio08"
## [19] "bio09"Resumir
Muitas vezes queremos fazer cálculos para todos as células de um raster. Podemos resumir informações de todos os pixels fazendo cálculos simples com todos os pixels de cada camada com a função raster::cellStats(), sendo x o argumento do objeto raster e stat o nome da função resumo, como “mean” ou “sum”.
# media de todas as celulas de altitude
raster::cellStats(x = ra_rc, stat = mean)## [1] 625.8273# media de todas as celulas de cada camada bioclimatica
raster::cellStats(x = st, stat = mean)## bio01 bio10 bio11 bio12 bio13 bio14
## -4.0378283 7.2035545 -13.8963286 550.0569022 93.4633916 15.3689993
## bio15 bio16 bio17 bio18 bio19 bio02
## 74.7084151 241.6525005 55.4149542 156.4237816 108.8950626 9.9432120
## bio03 bio04 bio05 bio06 bio07 bio08
## 34.5221528 880.1215546 13.9386423 -19.7938943 33.7325366 -0.9226276
## bio09
## -5.3774489Ou ainda, podemos analisar a frequência com que cada valor dos pixels ocorre, utilizando a função raster::freq().
# frequencia das celulas
raster::freq(x = ra_rc) %>% head()## value count
## [1,] 491 1
## [2,] 492 4
## [3,] 493 9
## [4,] 494 19
## [5,] 495 32
## [6,] 496 44# frequencia das celulas
raster::freq(x = st[[1]]) %>% head()## value count
## [1,] -55 319
## [2,] -54 4529
## [3,] -53 5778
## [4,] -52 6128
## [5,] -51 6090
## [6,] -50 7892Operações espaciais
As operações espaciais são modificações de objetos espaciais baseado em informações espaciais, como localização e formato. Seria impossível abordar todas as operações realizáveis, então demonstraremos as principais para dados vetoriais e raster.
Vetor
As principais operações espaciais para dados vetoriais são: 1) filtro espacial, 2) junção espacial, 3) agregação espacial e 4) distância espacial. Apresentaremos essas operações utilizando as princiais funções utilizando os dados de nascentes, hidrologia e cobertura da terra para o município de Rio Claro/SP.
Filtro espacial
Filtros espaciais são operações que realizam seleção de feições espaciais entre dois objetos espaciais (x e y). Existe uma grande quantidade de funções para realizar filtros espaciais no R, como podemos ver na Tabela (@ref(tab:tab-filtro-espacial), Pebesma & Bivand 2020). Essas funções verificam se cada feição em x mantém sua relação em y. Ao especificar o parâmetro sparse = FALSE, as funções retornam uma matriz lógica (comporta por TRUE e FALSE).
| Função | Descrição | Função inversa |
|---|---|---|
sf::st_contains() |
Nenhum dos pontos de x está fora de y | st_within |
sf::st_contains_properly() |
x contém y, e y não tem pontos em comum com a fronteira de x | NA |
sf::st_covers() |
Nenhum ponto de y se encontra no exterior de x | st_covered_by |
sf::st_covered_by() |
Inverso de sf::st_covers() |
NA |
sf::st_crosses() |
x e y têm alguns, mas não todos os pontos internos em comum | NA |
sf::st_disjoint() |
x e y não têm pontos em comum | st_intersects |
sf::st_equals() |
x e y são geometricamente iguais; o número de pedido dos nós pode ser diferente; idêntico a x contém y xND x dentro de y | NA |
sf::st_equals_exact() |
x e y são geometricamente iguais e têm ordem de nó idêntica | NA |
sf::st_intersects() |
x e y não são separados | st_disjoint |
sf::st_is_within_distance() |
x está mais perto de y do que uma determinada distância | NA |
sf::st_within() |
Nenhum dos pontos de y está fora de x | st_contains |
sf::st_touches() |
x e y têm pelo menos um ponto limite em comum, mas nenhum ponto interno | NA |
sf::st_overlaps() |
x e y têm alguns pontos em comum; a dimensão destes é idêntica à de x e y | NA |
sf::st_relate() |
Dado um padrão, retorna se x e y aderem a este padrão | NA |
Em nosso exemplo, utilizaremos a função sf::intersects() para filtrar as nascentes dentro de floresta para Rio Claro/SP. Essa função vai retornar a resposta binária se as nascentes estão (1) ou não (empty) dentro dos polígonos de floresta.
# filtro espacial
sf::st_intersects(x = rc_nas, y = rc_cob_floresta)## Sparse geometry binary predicate list of length 1220, where the predicate was `intersects'
## first 10 elements:
## 1: 1
## 2: 1
## 3: (empty)
## 4: 1
## 5: (empty)
## 6: (empty)
## 7: (empty)
## 8: (empty)
## 9: 1
## 10: (empty)Podemos usar essa mesma função em conjunto com a função dplyr::filter() para filtrar as nascentes dentro de florestas, mas agora com o argumento sparse = FALSE para valores lógicos funcionarem com o filtro.
# filtro espacial - interno
rc_nas_floresta_int <- rc_nas %>%
dplyr::filter(sf::st_intersects(x = ., y = rc_cob_floresta, sparse = FALSE))Ou ainda podemos utilizar o operador [] para realizar esse filtro, como podemos notar na Figura @ref(fig:fig-vetor-filtro-espacial-interno).
# filtro espacial com [] - interno
rc_nas_floresta_int <- rc_nas[rc_cob_floresta, ]plot(rc_2019_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(rc_cob_floresta$geometry, col = "forestgreen", add = TRUE)
plot(rc_nas_floresta_int$geometry, col = "blue", pch = 20, cex = 1, add = TRUE)
Nascentes dentro de florestas no município de Rio Claro/SP.
Entretanto, muitas vezes queremos fazer o filtro de feições que estão fora de feições de outro objeto espacial. Para isso, podemos usar a função sf::st_disjoint() ou ainda utilizando o operador [], mas com o argumento op, nesse caso utilizando a mesma função sf::st_disjoint() como operação (Figura @ref(fig:fig-vetor-filtro-espacial-externo)). Notar o segundo argumento vazio nesse filtro.
# filtro espacial - externo
rc_nas_floresta_ext <- rc_nas %>%
dplyr::filter(sf::st_disjoint(x = ., y = rc_cob_floresta, sparse = FALSE))
# filtro espacial com [] - externo
rc_nas_floresta_ext <- rc_nas[rc_cob_floresta, , op = st_disjoint]plot(rc_2019_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(rc_cob_floresta$geometry, col = "forestgreen", add = TRUE)
plot(rc_nas_floresta_ext$geometry, col = "steelblue", pch = 20, cex = 1, add = TRUE)
Nascentes fora de florestas no município de Rio Claro/SP.
Junção espacial
Outra operação muito usada dentro de análises geográficas é a junção espacial ou do inglês spatial join. A ideia base é muito semelhante com a junção baseada em atributos, mas aqui iremos atribuir o valor da tabela de atributos das feições de um objeto espacial y às feições que fazem intersecção com um objeto espacial x, de modo que esses valores sejam armazenados na tabela de atributos do primeiro objeto espacial.
Para exemplificar, vamos atribuir os valores dos polígonos de cobertura da terra aos pontos de nascentes para Rio Claro/SP, fazendo um agrupamento pela tabela de atributos para permitir criar o mapa da Figura @ref(fig:fig-vetor-juncao-espacial).
# juncao espacial
rc_nas_cob_jun <- rc_nas %>%
sf::st_join(x = ., y = rc_cob) %>%
dplyr::group_by(CLASSE_USO) %>%
dplyr::summarise(n = n())# plot
plot(rc_2019_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(rc_nas_cob_jun[1], col = c("blue", "orange", "gray30", "forestgreen", "green"), pch = 20, add = TRUE)
Junção espacial da cobertura da terra para as nascentes no município de Rio Claro/SP.
Agregação espacial
Muitas vezes queremos contabilizar quantas feições ou agregar valores de feições para polígonos. Podemos realizar essa operação usando as funções dplyr::group_by() e dplyr::summarise, ou utilizar a função aggregate(). Nesse exemplo, vamos contabilizar quantas nascentes há por cada polígono de cobertura da terra para o município de Rio Claro/SP (Figura @ref(fig:fig-vetor-agregacao-espacial)).
# agregacao espacial
rc_cob_nas_agre <- rc_nas %>%
aggregate(x = ., by = rc_cob, FUN = length)plot(rc_cob_nas_agre[1], axes = TRUE, graticule = TRUE, main = NA)
Agregação espacial contabilizando o número de nascentes para cada classe de cobertura da terra no município de Rio Claro/SP.
Distância espacial
A distância espacial é a distância calculada em duas dimensões (2D) entre um objeto espacial x e y baseado no CRS e para cada feição dos objetos espaciais. Para realizar esse cálculo, utilizamos a função sf::st_distance(). Em nosso exemplo, vamos calcular a distância das nascentes até a floresta mais próxima, e adicionando essa informação para cada ponto na tabela de atributos com a função dplyr::mutate(), para o município de Rio Claro/SP (Figura @ref(fig:fig-vetor-distancia-espacial)).
rc_nas_dist_flo <- rc_nas %>%
dplyr::mutate(dist_flo = sf::st_distance(rc_nas, rc_cob_floresta))plot(rc_2019_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(rc_cob_floresta$geometry, col = "forestgreen", add = TRUE)
plot(rc_nas_dist_flo[7], pch = 20, add = TRUE)
Distância espacial das nascentes até o fragmento de floresta mais próxima no município de Rio Claro/SP.
Raster
As principais operações espaciais para dados raster podem ser classificas, segundo @lovelace-etal-2019, em: 1) operações locais (por célula), 2) operações focais (por bloco de multiplas células regulares - e.g. 3x3), 3) operações zonais (por bloco de multiplas células irregulares) e 4) operações globais (por um ou vários rasters inteiros). Cada uma delas é aplicada para objetivos e escalas espaciais específicas. Para os exemplos desta seção, utilizaremos o dado raster de elevação para o município de Rio Claro/SP.
Operações locais
As operações locais contemplam todas as operações realizadas célula a célula em uma ou várias camadas de um objeto raster. A álgebra de raster é uma das mais comuns, simples e poderosas operações no R envolvendo rasters. Com ela podemos fazer operações simples através de operadores aritméticos (soma, subtração, multiplicação, divisão ou potenciação) entre dois ou mais objetos raster, ou utilizar funções para alterar todos os valores dos pixels como, por exemplo, as funções lo10() ou sqrt(), ou ainda a função raster::scale() para padronizar ou centralizar os valores dos rasters (Figura @ref(fig:fig-raster-local-aritmetico)).
# soma
ra_rc2 <- ra_rc + ra_rc
# log10
ra_rc_log10 <- log10(ra_rc)par(mfrow = c(1, 2))
raster::plot(ra_rc2, col = viridis::viridis(10))
plot(rc_2019$geom, col = NA, border = "red", lwd = 2, add = TRUE)
raster::plot(ra_rc_log10, col = viridis::viridis(10))
plot(rc_2019$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Rasters de soma e log10 do mapa de elevação para Rio Claro/SP.
par(mfrow = c(1, 1))Além das operação aritméticas, a álgebra de rasters também permite operações lógicas, como criar um novo raster (binário - composto por 1 quando a operação lógica é verdadeira, e 0 quanto é falsa). Em nosso caso, buscamos todos os pixels acima de 600 metros para o raster de elevação de Rio Claro/SP (Figura @ref(fig:fig-raster-local-logico)).
# acima de 600
ra_rc_acima_600 <- ra_rc > 600raster::plot(ra_rc_acima_600, col = viridis::viridis(10))
plot(rc_2019$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Operação local lógica mostrando todos os pixels acima de 600 metros de elevação para Rio Claro/SP.
Além dos operadores aritméticos, também podemos usar as funções raster::calc() (uma camada) e raster::overlay() (duas ou mais camadas) para realizar operações em todas as células. Elas funcionam com a criação de uma função específica através da função function(), para que esta seja aplicada em todas as células do raster. Essas funções são muito eficientes, portanto, são preferíveis para grandes conjuntos de dados raster. Exemplificaremos essa operação calculando o produto de todos os pixels por eles mesmos do raster de elevação de Rio Claro/SP (Figura @ref(fig:fig-raster-local-calc)).
# produto dos pixel - calc
ra_rc_prod <- raster::calc(x = ra_rc, fun = function(x){x * x})
ra_rc_prod ## class : RasterLayer
## dimensions : 370, 364, 134680 (nrow, ncol, ncell)
## resolution : 0.0008333333, 0.0008333333 (x, y)
## extent : -47.765, -47.46167, -22.55167, -22.24333 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 241081, 970225 (min, max)raster::plot(ra_rc_prod, col = viridis::viridis(10))
plot(rc_2019$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Operação local de multiplicação de todos os pixels por eles mesmos do raster de elevação para Rio Claro/SP.
A predição de objetos raster é outra aplicação extremamente útil de operações locais. A partir da relação entre variáveis respostas (e.g, pontos no espaço, como ocorrência ou riqueza de espécies), e variáveis preditoras (rasters contínuos de elevação, pH, precipitação, temperatura, cobertura da terra ou classe de solo), criamos modelos usando funções como lm(), glm(), gam() ou uma técnica de aprendizado de máquina, e fazemos predições espaciais aplicando os coeficientes estimados aos valores dos raster preditores (consulte a seção xx).
Por fim, a reclassificação de rasters é outra operação muito comum quando trabalhamos com rasters. Nela é realizada a classificação de intervalos de valores numéricos em grupos, e.g. agrupar um modelo digital de elevação em classes de valores. A função que faz essa operação é a raster::reclassify(). Ela possui dois argumentos: x que é o raster a ser reclassificado, e o segundo rcl para o qual devemos construir uma matriz de reclassificação, onde a primeira coluna é a extremidade inferior, a segunda coluna é a extremidade superior, e a terceira coluna representa o novo valor para os intervalos das colunas um e dois. Vamos reclassificar o raster de elevação de Rio Claro/SP para os intervalos 400–600, 600–800 e 800–1000 que são reclassificados para os valores 1, 2 e 3, respectivamente (Figura @ref(fig:fig-raster-local-reclassificacao)).
# matriz de reclassificacao
rcl <- matrix(c(400,600,1,
600,800,2,
800,1000,3),
ncol = 3, byrow = TRUE)
# reclassificao
ra_rc_rcl <- raster::reclassify(x = ra_rc, rcl = rcl)raster::plot(ra_rc_rcl, col = viridis::viridis(3))
plot(rc_2019$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Operação local de reclassificação para três classes de elevação para Rio Claro/SP.
Operações focais
As operações focais levam em consideração uma célula central e seus vizinhos. A vizinhança (também chamada de janela móvel - moving window) tipicamente é composta de células de 3 por 3 (célula central e seus oito vizinhos), mas pode assumir outra forma. A operação focal aplica uma função de agregação a todas as células dentro da vizinhança especificada, e usa a saída correspondente como o novo valor para a célula central, e segue para a próxima célula central e seus vizinhos. Essa operação é realizada através da função raster::focal(). O parâmetro x especifica o raster de entrada, o parâmetro w define a janela móvel por uma matriz cujos valores correspondem a pesos, e por fim o parâmetro fun especifica a função que desejamos aplicar às céluas, como min(), max(), sum(), mean(), sd() ou var(). Existem diversas aplicações dessa operação para dados raster, como no processamento de imagens de satélite (ver mais em Wegmann et al. 2016). Outra utilizade é para o cálculo de características topográficas, como declividade, aspecto e direções de fluxo. Para calcular essas métricas específicas, podemos utilizar a função raster::terrain().
Para nosso exemplo, vamos realizar o cálculo do desvio padrão da elevação e a métrica de aspecto (orientação da vertente) para o raster de elevação em Rio Claro/SP (Figura @ref(fig:fig-raster-focal)).
# janela movel
ra_rc_focal_sd <- raster::focal(x = ra_rc, w = matrix(data = 1, nrow = 3, ncol = 3), fun = sd)
# declividade
ra_rc_asp <- raster::terrain(x = ra_rc, opt = "aspect")par(mfrow = c(1, 2))
raster::plot(ra_rc_focal_sd, col = viridis::viridis(10))
plot(rc_2019$geom, col = NA, border = "red", lwd = 2, add = TRUE)
raster::plot(ra_rc_asp, col = viridis::viridis(10))
plot(rc_2019$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Cálculo do desvio padrão da elevação para uma janela de 3x3 e do aspecto para Rio Claro/SP.
par(mfrow = c(1, 1))Operações zonais
As operações zonais aplicam uma função de agregação a várias células de um raster. Geralmente usa-se um segundo raster categórico para definir as zonas, de modo que as células raster que definem a zona não precisam ser vizinhas, como na operação focal. O resultado de uma operação zonal é uma tabela de resumo agrupada por zona, explicando porque essa operação também é conhecida como estatística zonal. Isso é um contraste com as operações focais que retornam um objeto raster.
A operação zonal é realizada através da função raster::zonal(), que recebe de entrada no argumento x o raster contínuo, em z o raster categórico, e em fun a função que irá resumir as células. Em nosso exemplo, vamos calcular diversas medidas resumo da elevação com a função summary() para cada classe de elevação que criamos anteriormente.
# estatistica zonal
ra_rc_zonal <- data.frame(raster::zonal(ra_rc, ra_rc_rcl, fun = "summary"))
colnames(ra_rc_zonal) <- c("zona", "min", "1qt", "mediana", "media", "3qt", "max")
ra_rc_zonal## zona min 1qt mediana media 3qt max
## 1 1 491 552 574.0 567.5995 589 600
## 2 2 601 620 640.0 650.6829 670 800
## 3 3 801 817 832.5 834.2732 846 985Operações globais
As operações globais usam todo o conjunto de dados raster representando uma única zona. As operações globais mais comuns são estatísticas descritivas para todos os pixels do raster, utilizando a função raster::cellStats() ou raster::freq(), já vistas. Além das estatísticas descritivas, podemos gerar rasters de distância, que calcula a distância de cada célula a uma ou um grupo células-alvo específica, utilizando a função raster::distance().
Em nosso exemplo, vamos selecionar os pixels abaixo de 500 m do raster de elevação e calcular a distância Euclidiana (Figura @ref(fig:fig-raster-global)).
# distancia
ra_rc_abaixo_500 <- raster::calc(x = ra_rc, fun = function(x) ifelse(x < 500, 1, NA))
ra_rc_global_dist <- raster::distance(ra_rc_abaixo_500)plot(ra_rc_global_dist, col = viridis::viridis(10))
plot(ra_rc_abaixo_500, add = TRUE, col = "white", legend = FALSE)
plot(rc_2019$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Raster de distância Euclidiana dos pixels abaixo de 500 m de elevação para Rio Claro/SP.
Operações geométricas
As operações geométricas realizam modificações em objetos espaciais baseado na geometria do vetor ou do raster e na interação e conversão entre vetor-raster.
As operações geométricas vetoriais podem ser unárias funcionam em uma única geometria, ou binárias que modificam uma geometria com base na forma de outra. Ainda podemos fazer transformações para alterar os tipos vetores, que irá refletir se as feições são únicas ou múltiplas, inclusive na tabela de atributos.
As operações geométricas em rasters envolvem mudar a posição, tamanho e número dos pixels subjacentes e atribuir-lhes novos valores. Por fim, podemos ainda fazer operações de interações e conversões entre raster-vetor para ajustar rasters à vetores, assim como converter um objeto espacial vetorial para raster e vice-versa.
Vetor
Como dissemos, as operações geométricas em vetores irão criar ou alterar a geometria de objetos da classe sf, podendo fazer alterações em única geometria (unárias): 1) simplificação, 2) centróides, 3) pontos aleatórios, 4) buffers, 5) polígono convexo, 7) polígonos de Voronoi, 7) quadrículas e hexágonos; ou modificar uma geometria com base na forma de outra geometria (binárias): 8) união e 9) recortes; ou ainda fazer transformações de tipo.
Para exemplificar as operações geométricas com vetores, vamos utilizar os dados do limite, nascentes, hidrologia e cobertura da terra para o município de Rio Claro/SP.
Simplificação
A simplificação possui o intuito de generalizar linhas ou polígonos, diminuindo assim suas complexidades em relação ao número de vértices. É utilizada para representação em mapas menores ou mapas interativos (seção xx), ou ainda quando um objeto vetorial é muito grande. A função utilizada é a sf::st_simplify(), que usa o argumento x para uma geometria de entrada e dTolerance para controlar o nível de generalização nas unidades do mapa. Em nosso exemplo, simplificaremos a hidrografia de Rio Claro/SP (Figura @ref(fig:fig-vetor-simplificacao)).
# simplificacao
rc_hid_simplificado <- sf::st_simplify(x = rc_hid, dTolerance = 1000)plot(rc_2019_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(rc_hid$geometry, col = "steelblue", lwd = 2, add = TRUE)
plot(rc_hid_simplificado$geometry, col = adjustcolor("black", .7), add = TRUE)
Simplificação da hidrografia para Rio Claro/SP.
Centroides
A operação de centroides identifica o centro de objetos geográficos, geralmente o centro de massa das feições. É utilizado para gerar um ponto simples para representações complexas ou para estimar a distância entre polígonos utilizando esse centroide. Podemos calculá-los com a função sf::st_centroids(), ou com a função sf::st_point_on_surface() para garantir que esses centroides caiam dentro das geometrias. Aqui calcularemos o centroide do município de Rio Claro/SP (Figura @ref(fig:fig-vetor-centroide)).
# centroides
rc_2019_sirgas2000_utm23s_cent <- sf::st_centroid(rc_2019_sirgas2000_utm23s)## Warning in st_centroid.sf(rc_2019_sirgas2000_utm23s): st_centroid assumes
## attributes are constant over geometries of xplot(rc_2019_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(rc_2019_sirgas2000_utm23s_cent$geom, cex = 3, pch = 20, add = TRUE)
Centroide do limite do município de Rio Claro/SP.
Pontos aleatórios
Por vezes precisamos criar algum padrão aleatoriotório dentro de um contexto espacial. Isso pode ser realizado de diversas formas. Uma delas é a criação de pontos aleatórios dentro de um polígono. Podemos realizar essa operação com a função sf::st_sample(). Para essa função, dois argumentos são utilizados: x uma geometria de entrada e o size indicando o número de pontos à ser criado. Outro argumento bastante interessante é o type, indicando o tipo de amostragem espacial (aleatório, regular ou triangular). Para nosso exemplo, vamos fixar a amostragem utilizando a função set.seed() e sortear 30 pontos para o limite do município de Rio Claro/SP.
# fixar amostragem
set.seed(42)
# pontos aleatoriotorios
rc_2019_sirgas2000_utm23s_pontos_aleatorios <- sf::st_sample(rc_2019_sirgas2000_utm23s, size = 30)plot(rc_2019_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(rc_2019_sirgas2000_utm23s_pontos_aleatorios, pch = 20, add = TRUE)
Sorteio de 30 pontos aleatoriotório para Rio Claro/SP.
Buffer
Buffers são polígonos que representam a área dentro de uma determinada distância de um elemento geométrico, independentemente de ser um ponto, linha ou polígono. O buffer é comumente utilizado para análise de dados geográficos, geralmente sendo entendio como uma unidade amostral, delimitando uma porção no entorno de algum elemento ou evento, como as condições climáticas ou da estrutura da paisagem para uma amostragem, ou as características de cobertura da terra ao longo de um corpo d’água, e.g. a Área de Preservação Permanente (APP).
A função utilizada para criar buffers é a sf::st_buffer(), que requer pelo menos dois argumentos: x uma geometria de entrada e o dist uma distância para o buffer, fornecido nas unidades do CRS da geometria de entrada. Em nosso exemplo, vamos criar buffers de 1000 metros para os 30 pontos aleatórios criados anteriormente para o município de Rio Claro/SP (Figura @ref(fig:fig-vetor-buffer)).
# buffer
rc_2019_sirgas2000_utm23s_pontos_aleatorios_buffer <- sf::st_buffer(x = rc_2019_sirgas2000_utm23s_pontos_aleatorios, dist = 1000)plot(rc_2019_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(rc_2019_sirgas2000_utm23s_pontos_aleatorios_buffer, col = NA, lwd = 2, border = "red", add = TRUE)
plot(rc_2019_sirgas2000_utm23s_pontos_aleatorios, pch = 20, cex = 1, add = TRUE)
Buffers de 1000 metros para os 30 pontos aleatórios no município de Rio Claro/SP.
Polígono convexo
Uma análise bastante comum, principalmente realizada pela IUCN, é a criação de polígonos convexos, para definir a extensão de ocorrência de uma espécie (Extent of occurrence - EOO). Nesse sentido, essa operação irá ligar os pontos externos de um conjunto de pontos e criar um polígono a partir deles. Podemos criar esse polígono com a função sf::st_convex_hull(). Um único passo que precisamos adiantar é utilizar a função sf::st_union() para unir todos os pontos e criar um objeto sf MULTIPOINT, já iremos explicado com mais detalhes adiante. Vamos utilizar os pontos aleatórios que criamos anteriormente para criar o polígono convexo (Figura @ref(fig:fig-vetor-convexo)).
# poligono convexo
rc_2019_sirgas2000_utm23s_convexo <- rc_2019_sirgas2000_utm23s_pontos_aleatorios %>%
sf::st_union() %>%
sf::st_convex_hull()plot(rc_2019_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(rc_2019_sirgas2000_utm23s_convexo, col = NA, lwd = 2, border = "red", add = TRUE)
plot(rc_2019_sirgas2000_utm23s_pontos_aleatorios, pch = 20, cex = 1, add = TRUE)
Polígono convexo para os 30 pontos criados aleatoriotoriamente para Rio Claro/SP.
Polígonos de Voronoi
Uma outra forma de criar polígonos para resumir dados espaciais é através de Polígonos de Voronoi ou Diagrama de Voronoi. Nele, polígonos irregulares são criados à partir da proximidade de pontos, de modo a estimar uma área de abrangência no entorno dos mesmos (Okabe et al. 2000). Esses polígonos podem ser criados com a função sf::st_voronoi(), mas precisamos novamente utilizar a função sf::st_union() para unir todos os pontos e criar um objeto sf MULTIPOINT. Vamos utilizar os pontos aleatórios que criamos anteriormente para criar o polígono de Voronoi (Figura @ref(fig:fig-vetor-voronoi)).
# poligono de voronoi
rc_2019_sirgas2000_utm23s_voronoi <- rc_2019_sirgas2000_utm23s_pontos_aleatorios %>%
sf::st_union() %>%
sf::st_voronoi()plot(rc_2019_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(rc_2019_sirgas2000_utm23s_voronoi, col = NA, lwd = 2, border = "red", add = TRUE)
plot(rc_2019_sirgas2000_utm23s_pontos_aleatorios, pch = 20, cex = 1, add = TRUE)
Polígonos de Voronoi para os 30 pontos criados aleatoriotoriamente para Rio Claro/SP.
Quadrículas e hexágonos
Muitas vezes precisamos criar unidades espaciais idênticas e igualmente espaçadas para resumir informações dispersas por toda a nossa área de estudo. Uma prática muito comum é a criação de um gride de pontos ou quadrículas em toda a área de estudo, e depois utilizar essas geometrias para associar ou resumir informações espacializadas (exemplo na seção xx), como a IUCN utiliza para a análise de área de ocupação (Are of occupancy - AOO). Além das quadrículas, uma outra geometria que se tornou bastante comum para as finalidades descritas, é a criação de hexágonos, que além de serem mais esteticamente atraentes, possuem uma explicação matemática de sua melhor funcionadade para análises espaciais em Ecologia (Birch et al. 2007).
A função utilizada para criar esses grides é a sf::st_make_grid(), que requer pelo menos dois argumentos: x uma geometria de entrada e o cellsize indicando o tamanho do gride a ser criado, fornecido nas unidades do CRS da geometria de entrada. Há diversos outros argumentos, mas os mais importantes são o square que irá definir se o gride será de quadriculas ou de hexágonos, e o what que irá definir se iremos gerar polígonos, cantos ou centroides.
Em nosso exemplo, vamos criar quadrículas e hexágonos de 2000 metros (i.e. 4000000 metros quadrados) para o município de Rio Claro/SP (Figura @ref(fig:fig-vetor-quad) e Figura @ref(fig:fig-vetor-hex)). Podemos ainda utilizar as funções de filtros espaciais (Tabela @ref(tab:tab-filtro-espacial)) para definir como iremos selecionar esses elementos para a área de estudo. Aqui utilizamos a função sf::st_intersects().
# quadriculas
rc_2019_sirgas2000_utm23s_grid <- sf::st_make_grid(x = rc_2019_sirgas2000_utm23s, cellsize = 2000, what = "polygons") %>%
sf::st_as_sf() %>%
dplyr::filter(sf::st_intersects(x = ., y = rc_2019_sirgas2000_utm23s, sparse = FALSE))
# centroides das quadriculas
rc_2019_sirgas2000_utm23s_grid_cent <- rc_2019_sirgas2000_utm23s %>%
sf::st_make_grid(cellsize = 2000, what = "centers") %>%
sf::st_as_sf() %>%
dplyr::filter(sf::st_intersects(x = ., y = sf::st_union(rc_2019_sirgas2000_utm23s_grid), sparse = FALSE))plot(rc_2019_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(rc_2019_sirgas2000_utm23s_grid, col = NA, border = "red", lwd = 2, add = TRUE)
plot(rc_2019_sirgas2000_utm23s_grid_cent, pch = 20, add = TRUE)
Quadrículas de 2000 metros e centroides para Rio Claro/SP.
# hexagonos
rc_2019_sirgas2000_utm23s_hex <- rc_2019_sirgas2000_utm23s %>%
sf::st_make_grid(cellsize = 2000, square = FALSE) %>%
sf::st_as_sf() %>%
dplyr::filter(sf::st_intersects(x = ., y = rc_2019_sirgas2000_utm23s, sparse = FALSE))
# centroides de hexagonos
rc_2019_sirgas2000_utm23s_hex_cent <- rc_2019_sirgas2000_utm23s %>%
sf::st_make_grid(cellsize = 2000, square = FALSE, what = "centers") %>%
sf::st_as_sf() %>%
dplyr::filter(sf::st_intersects(x = ., y = sf::st_union(rc_2019_sirgas2000_utm23s_hex), sparse = FALSE))plot(rc_2019_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(rc_2019_sirgas2000_utm23s_hex, col = NA, border = "red", lwd = 2, add = TRUE)
plot(rc_2019_sirgas2000_utm23s_hex_cent, pch = 20, add = TRUE)
Hexágonos de 2000 metros e centroides para Rio Claro/SP.
União (“dissolver”)
Como vimos na seção xx, a agregação por atributos podemos dissolver as geometrias de polígonos do mesmo grupo pelos valores da tabela de atributos, onde, naquele exemplo, contabilizamos quantas nascentes havia por polígono de cobertura da terra para o município de Rio Claro/SP (Figura @ref(fig:fig-vetor-agregacao-espacial)).
Nesta seção, vamos utilizar a função sf::st_union() para unir diversas feições em uma só, dissolvendo os limites entre elas. Vamos utilizar de exemplo os buffers que criamos a partir dos 30 pontos aleatórios (Figura @ref(fig:fig-vetor-uniao)).
# uniao
rc_2019_sirgas2000_utm23s_pontos_aleatorios_buffer_uniao <- sf::st_union(rc_2019_sirgas2000_utm23s_pontos_aleatorios_buffer)plot(rc_2019_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(rc_2019_sirgas2000_utm23s_pontos_aleatorios_buffer_uniao, col = adjustcolor("blue", .1), add = TRUE)
União - dissolução - dos buffers criados a partir de pontos aleatórios para Rio Claro/SP.
Recorte (“clipar”)
O recorte realiza um subconjunto espacial envolvendo dois objetos espaciais. O recorte é aplicado somente a linhas e polígonos, ou seja, usaremos linhas e polígonos para recortar linhas ou polígonos. Esse recorte pode ser realizado de três formas: 1) intersecção (subconjunto das geometrias sobrepostas entre os dois objetos), 2) diferença (subconjunto das geometrias do primeiro objeto sem sobreposição com o segundo objeto), e 3) diferença simétrica (apenas as geometrias não sobrepostas entre os dois objetos). Respectivamente para cada uma dessas operações temos funções específicas: sf::st_intersection(), sf::st_difference() e sf::st_sym_difference().
Para nosso exemplo, faremos o recorte da hidrogradia em relação aos buffers criados e unidos para os 30 pontos aleatórios em Rio Claro/SP. Primeiramente, iremos fazer o recorte para dentro dos buffers com a função sf::st_intersection() (Figura @ref(fig:fig-vetor-interseccao)).
# recorte - interseccao
rc_hid_interseccao <- sf::st_intersection(x = rc_hid, y = rc_2019_sirgas2000_utm23s_pontos_aleatorios_buffer_uniao)## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriesplot(rc_2019_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(rc_2019_sirgas2000_utm23s_pontos_aleatorios_buffer_uniao, col = adjustcolor("blue", .1), add = TRUE)
plot(rc_hid_interseccao$geometry, col = "blue", add = TRUE)
Recorte da hidrografia para dentro dos buffers dos 30 aleatórios para Rio Claro/SP.
Para nosso segundo exemplo, realizamos o recorte da hidrogradia em relação aos buffers, mas agora para fora dos buffers utilizando a função sf::st_difference() (Figura @ref(fig:fig-vetor-interseccao)).
# recorte - diferenca
rc_hid_diferenca <- sf::st_difference(x = rc_hid, y = rc_2019_sirgas2000_utm23s_pontos_aleatorios_buffer_uniao)## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriesplot(rc_2019_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(rc_2019_sirgas2000_utm23s_pontos_aleatorios_buffer_uniao, col = adjustcolor("blue", .1), add = TRUE)
plot(rc_hid_diferenca$geometry, col = "blue", add = TRUE)
Recorte da hidrografia para fora dos buffers dos 30 aleatórios para Rio Claro/SP.
Transformações de tipo
Esse tópico possui muitas funcionalidades, que são exploradas no tópico “5.2.7 Type transformations” de @lovelace-etal-2019. Aqui, nosso interesse principal é em relação à transformação dos tipos de objetos espaciais da classe sf: MULTIPOINT, MULTILINESTRING e MULTIPOLYGON, para POINT, LINESTRING e POLYGON. Muitas vezes as feições de nossos objetos, i.e., as linhas da tabela de atributos, estão agrupadas em apenas um linha da tabela. Quando o objeto espacial está nesse formato, geralmente em alguma classe dessas (MULTIPOINT, MULTILINESTRING e MULTIPOLYGON), não temos como realizar operações espaciais ou geométricas para cada feição, e precisamos separá-las cada uma em uma linha para que operações como o cálculo de comprimento ou área seja possível para cada feição (veja seção xx).
Dessa forma, podemos utilizar a função sf::st_cast() para fazer essas transformações e atribuir cada feição à uma linha da tabela de atributos. Como exemplo, vamos separar os fragmentos de floresta e calcular a área para cada feição em hectares (Figura @ref(fig:fig-vetor-tipo)).
# transformacao de tipo
rc_cob_floresta_polygon <- rc_cob_floresta %>%
sf::st_cast("POLYGON") %>%
dplyr::mutate(area_ha = sf::st_area(.)/1e4 %>% round(2))## Warning in st_cast.sf(., "POLYGON"): repeating attributes for all sub-geometries
## for which they may not be constantplot(rc_2019_sirgas2000_utm23s$geom, col = "gray", main = NA, axes = TRUE, graticule = TRUE)
plot(rc_cob_floresta_polygon["area_ha"], col = viridis::viridis(100), add = TRUE)## Warning in plot.sf(rc_cob_floresta_polygon["area_ha"], col =
## viridis::viridis(100), : col is not of length 1 or nrow(x): colors will be
## recycled; use pal to specify a color palette
Transformação do vetor de florestas em POLYGON e cálculo da área para cada feição para Rio Claro/SP.
Raster
As operações geométricas em rasters envolvem mudar a posição, tamanho e número dos pixels e atribuir novos valores, geralmente aumentando ou diminuindo o tamanho desses pixels. Essas operações permitem alinhar rasters de diversas fontes, fazendo com que compartilhem uma correspondência entre os pixels, permitindo que eles sejam processados todos juntos, ou simplesmente permita a realização de análises que demorariam muito, caso os rasters possuam um tamanho de pixel muito pequeno. Importante frisar que essas operação funcionam para as três classes dos objetos raster: RasterLayer, RasterBrick e RasterStack.
Para exemplificar as operações geométricas com rasters, vamos utilizar os dados de elevação para o município de Rio Claro/SP e bioclimáticos para o mundo.
Agregação
Na agregação de rasters iremos aumentar o tamanho dos pixels (diminuindo a resolução), agregando os valores dos pixels em um pixel maior. Podemos realizar essa operação com a função raster::aggregate(), que possui três argumentos: x corresponde ao objeto raster de entrada, fact é o fator de agregação e corresponde ao número que definirá o novo tamanho do pixel (e.g. se um raster tem resolução de 90 m, um fator de agregação de 10 fará com o novo raster tenha a resolução de 900 m), e fun é a função utilizada para realizar a agregação dos pixels (Figura @ref(fig:fig-raster-agregacao)).
Em nosso exemplo, vamos aumentar o tamanho dos pixels para 900 metros do raster de elevação para Rio Claro/SP.
# agregacao - aumentar o tamanho do pixel
ra_rc_sirgas2000_utm23s_agre_media <- raster::aggregate(x = ra_rc_sirgas2000_utm23s, fact = 10, fun = "mean")
ra_rc_sirgas2000_utm23s_agre_media## class : RasterLayer
## dimensions : 40, 37, 1480 (nrow, ncol, ncell)
## resolution : 900, 900 (x, y)
## extent : 214554.4, 247854.4, 7502625, 7538625 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=23 +south +ellps=GRS80 +units=m +no_defs
## source : memory
## names : srtm_27_17
## values : 506.0024, 922.8709 (min, max)raster::plot(ra_rc_sirgas2000_utm23s_agre_media, col = viridis::viridis(10))
plot(rc_2019_sirgas2000_utm23s$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Agregação (aumento do pixel para 900 metros) utilizando a média para o raster de elevação para Rio Claro/SP.
Desagregação
De modo contrátio, na desagregação de rasters iremos diminuir o tamanho dos pixels (aumentando a resolução), preenchendo com novos valores. Podemos realizar essa operação com a função raster::desaggregate(), que assim como a função anterior, possui três argumentos: x corresponde ao objeto raster de entrada, fact é o fator de desagregação e corresponde ao número que definirá o novo tamanho do pixel (e.g. se um raster tem resolução de 90 m, um fator de desagregação de 2 fará com que o novo raster tenha a resolução de 9 m), e method é a função utilizada para realizar a desagregação dos pixels (Figura @ref(fig:fig-raster-desagregacao)).
Nesso exemplo, vamos diminuir o tamanho dos pixels para 45 metros do raster de elevação para Rio Claro/SP.
# desagregacao - diminuir o tamanho do pixel
ra_rc_desg_bil <- raster::disaggregate(x = ra_rc_sirgas2000_utm23s, fact = 2, method = "bilinear")
ra_rc_desg_bil## class : RasterLayer
## dimensions : 792, 728, 576576 (nrow, ncol, ncell)
## resolution : 45, 45 (x, y)
## extent : 214554.4, 247314.4, 7502985, 7538625 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=23 +south +ellps=GRS80 +units=m +no_defs
## source : memory
## names : srtm_27_17
## values : 493.7046, 986.5187 (min, max)raster::plot(ra_rc_desg_bil, col = viridis::viridis(10))
plot(rc_2019_sirgas2000_utm23s$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Desagregação (diminuição do pixel para 45 metros) utilizando o método bilinear para o raster de elevação para Rio Claro/SP.
Alinhamento de rasters
Muitas vezes queremos ir além de ajustar o tamanho do pixel, ajustando também a extensão, número e origem dos pixels para várias camadas rasters, principalmente se precisamos criar objetos das classes RasterBrick ou RasterStack. Dessa forma, podemos utilizar a função raster::compareRaster() para comparar os rasters em relaçao a extensão, número de linhas e colunas, projeção, resolução e origem (ou um subconjunto dessas comparações).
Podemos utilizar a função raster::resample() para fazer esse alinhamento, ou ainda a função gdalUtils::align_rasters(). Para a primeira função, os argumentos são x para o raster de entrada, y para o raster de alinhamento e method para o método utilizado no alinhamento. Para nosso exemplo, vamos ajustar uma camada bioclimática (BIO01) à camada de elevação para Rio Claro/SP (Figura @ref(fig:fig-raster-reamostragem)).
# reamostragem
st_rc <- raster::resample(x = st$bio01, y = ra_rc, method = "bilinear")
st_rc## class : RasterLayer
## dimensions : 370, 364, 134680 (nrow, ncol, ncell)
## resolution : 0.0008333333, 0.0008333333 (x, y)
## extent : -47.765, -47.46167, -22.55167, -22.24333 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : bio01
## values : 19.8383, 20.60492 (min, max)raster::plot(st_rc, col = viridis::viridis(10))
plot(rc_2019_sirgas2000_utm23s$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Reamostragem (alinhamento dos rasters) utilizando o método bilinear para alinhar o raster bioclimático (BIO01) ao de elevação para Rio Claro/SP.
Interações raster-vetor
Podemos fazer operações da interação entre objetos vetoriais e raster, como ajustes da extensão e limite do raster para vetores (corte e máscara), extração dos valores dos pixels para vetores (pontos, linhas e polígonos), e estatísticas zonais dos valores dos pixels dos raster para um vetor (polígonos).
Cortes e máscaras
Muitas vezes precisamos ajustar o tamanho de um objeto raster à uma área menor de interesse, geralmente definido por um objeto vetorial. Para realizar essa operação, dispomos de duas funções: raster::crop() e raster::mask(), sendo que ambos os objetos raster a ser reduzido e vetor como molde precisam estar no mesmo CRS.
A função raster::crop() ajusta o raster à extensão do vertor. Como exemplo, vamos retomar o raster de elevação original baixado do site xx e importado na seção xx (Figura @ref(fig:fig-raster-dem)). Primeiramente, vamos usar a função raster::crop() para ajustar esse raster à extensão do limite do município de Rio Claro/SP (Figura @ref(fig:fig-raster-corte)).
# crop - adjuste da extensão
ra_rc_crop <- raster::crop(ra, rc_2019)raster::plot(ra_rc_crop, col = viridis::viridis(10))
plot(rc_2019$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Ajuste da extensão do raster de elevação para a extensão de Rio Claro/SP.
Para ajustar o raster ao limite do município de Rio Claro/SP, vamos usar a função raster::mask(). É importante notar que essa função preenche com NAs os pixels que estão fora do limite do polígono e não ajusta a extensão (Figura @ref(fig:fig-raster-corte)).
# mask - adjuste ao limite
ra_rc_mask <- raster::mask(ra, rc_2019)## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
## definition## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Sistema de Referencia Geocentrico para las
## AmericaS 2000 in Proj4 definitionraster::plot(ra_rc_mask, col = viridis::viridis(10))
plot(rc_2019$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Ajuste do raster de elevação para o limite de Rio Claro/SP.
Para ajustar o raster à extensão e ao limite do município de Rio Claro/SP, precisamos utilizar conjuntamente as funções raster::crop() e raster::mask() (Figura @ref(fig:fig-raster-corte-mascara)).
# crop e mask - ajuste da extensão e do limite
ra_rc_crop_mask <- ra %>%
raster::crop(rc_2019) %>%
raster::mask(rc_2019)## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
## definition## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Sistema de Referencia Geocentrico para las
## AmericaS 2000 in Proj4 definitionraster::plot(ra_rc_crop_mask, col = viridis::viridis(10))
plot(rc_2019$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Ajuste da extensão e do limite do raster de elevação para Rio Claro/SP.
A função raster::mask() possui ainda um argumento chamado inverse, que cria uma máscara inversa ao limite, preenchendo com NA o pixels internos ao limite do polígono, como podemos ver para o raster de elevação e o limite de Rio Claro/SP (Figura @ref(fig:fig-raster-corte-mascara-inverso)).
# crop e mask inversa - ajuste da extensão e do limite inverso
ra_rc_crop_mask_inv <- ra %>%
raster::crop(rc_2019) %>%
raster::mask(rc_2019, inverse = TRUE)## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
## definition## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Sistema de Referencia Geocentrico para las
## AmericaS 2000 in Proj4 definitionraster::plot(ra_rc_crop_mask_inv, col = viridis::viridis(10))
plot(rc_2019$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Ajuste da extensão e do limite externo do raster de elevação para Rio Claro/SP.
Extração
A interação entre raster-vetor de extração é o processo que identifica e retorna valores associados de pixels de um raster com base em um objeto vetorial. É uma operação extramamente comum em análises geográficas, principalmente para associar valores de raster ambientais (contínuos ou categóricos) a pontos de ocorrência ou amostragem. Os valores retornados irão depender do tipo vetor (pontos, linhas ou polígonos) e de argumentos da função raster::extract() que alteram o funcionamento da extração.
Em nosso exemplo, vamos extrair os valores do raster de elevação para as nascentes do município de Rio Claro/SP (Figura @ref(fig:fig-raster-extracao)).
# extracao
rc_nas_ele <- rc_nas %>%
dplyr::mutate(elev = raster::extract(x = ra_rc_sirgas2000_utm23s, y = .))## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
## definition## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Sistema de Referencia Geocentrico para las
## AmericaS 2000 in Proj4 definition# plot
plot(rc_nas_ele["elev"], pch = 20, main = NA, axes = TRUE, graticule = TRUE)
Extração dos valores de elevação para as nascentes de Rio Claro/SP.
Além da extração dos valores totais, podemos resumir os valores dos pixels com a mesma operação de extração, utilizando ainda a função raster::extract(), mas utilizando uma função para resumir os valores dos pixels para um polígono, operação também denominada de estatística zonal (agora para vetores). Já vimos que ela pode ser realizada entre rasters na seção xx, mas aqui a realizaremos para rasters e vetores.
Para o exemplo, vamos calcular a elevação média dos valores para os hexágonos que criamos para o limite de Rio Claro/SP( Figura @ref(fig:fig-raster-zonas)).
# extracao - estatistica por zonas
rc_2019_sirgas2000_utm23s_hex_alt <- rc_2019_sirgas2000_utm23s_hex %>%
dplyr::mutate(elev_mean = raster::extract(x = ra_rc_sirgas2000_utm23s, y = rc_2019_sirgas2000_utm23s_hex, fun = mean, na.rm = TRUE))## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
## definition## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Sistema de Referencia Geocentrico para las
## AmericaS 2000 in Proj4 definition# plot
plot(rc_2019_sirgas2000_utm23s_hex_alt["elev_mean"], pch = 20, main = NA, axes = TRUE, graticule = TRUE)
Extração dos valores de elevação e resumo pela média para os hexágonos de Rio Claro/SP.
Conversões raster-vetor
Por fim, podemos ainda fazer operações de conversão entre objetos vetoriais para raster e vice-versa. Nessas operações, podemos resumir ou transformar objetos vetoriais (pontos, linhas ou polígonos) para rasters, escolhendo um raster previamente existente, processo denominado rasterização. Também podemos realizar o processo inverso, i.e., transformar o raster em um vetor, podendo esse vetor ser um gride pontos, linhas ou polígonos, operação chamada de vetorização.
Rasterização
A conversão de vetor para raster pode ser realizada de pontos para rasters. Nesse processo, podemos utilizar uma função para resumir os dados pontuais para os pixels do raster que iremos criar. Para essa operação, podemos utilizar a função raster::rasterize(), com o argumento x sendo o vetor de pontos de entrada, y o raster base, field a coluna ou campo da tabela de atributos do ojeto vetorial para os quais os valores serão utilizados e fun a função utilizada para agregação dos dados.
Aqui, vamos contabilizar a quantidade de nascentes por pixel, utilizando como base o raster para o qual mudamos a resolução para 900 metros (@ref(fig:fig-raster-rasterizacao-pontos)).
# rasterizar pontos
rc_nas_rasterizacao <- raster::rasterize(x = rc_nas, y = ra_rc_sirgas2000_utm23s_agre_media, field = 1, fun = "count")## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
## definition## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Sistema de Referencia Geocentrico para las
## AmericaS 2000 in Proj4 definitionraster::plot(rc_nas_rasterizacao, col = viridis::viridis(10))
plot(rc_nas$geometry, pch = 20, cex = .5, col = adjustcolor("gray", .5), add = TRUE)
Rasterização das nascentes, com a operação de contabilização de pontos para Rio Claro/SP.
Além de pontos, podemos também rasterizar linhas. Aqui vamos contanilizar as linhas da hidrografia simplificada para Rio Claro/SP (Figura @ref(fig:fig-raster-rasterizacao-linhas)).
# rasterizar linhas
rc_hid_rasterizacao <- raster::rasterize(x = rc_hid_simplificado,
y = ra_rc_sirgas2000_utm23s_agre_media,
field = 1, fun = "count")## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
## definition## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Sistema de Referencia Geocentrico para las
## AmericaS 2000 in Proj4 definitionraster::plot(rc_hid_rasterizacao, col = viridis::viridis(10))
plot(rc_hid_simplificado$geom, col = "gray", add = TRUE)
Rasterização da hidrografia, com a operação de contabilização de linhas para Rio Claro/SP.
Podemos ainda rasterizar polígonos, de modo que cada pixel do raster a ser criado irá receber o valor da tabela de atributos, ou uma análise pelo vizinho mais próximo no caso de um campo categórico, como a cobertura da terra, que também vai depender da resolução do raster base e do tamanho da feição do polígono. Para nosso exemplo, antes de criar o raster vamos transforma a coluna de classe de cobertura da terra em factor (Figura @ref(fig:fig-raster-rasterizacao-poligonos)). Entretanto, essa operação de rasterização tente a demorar muito no caso de polígonos muito detalhados ou um raster com pixels muito pequenos, sendo que dois pacotes aceleram esse processamento (fasterize e gdalUtils), com suas respectivas funções: fasterize::fasterize() e gdalUtils::gdal_rasterize().
# rasterizar poligonos
rc_cob_rasterizacao <- rc_cob %>%
dplyr::mutate(classe = as.factor(CLASSE_USO)) %>%
raster::rasterize(x = ., y = ra_rc_sirgas2000_utm23s_agre_media, field = "classe")## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
## definition## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Sistema de Referencia Geocentrico para las
## AmericaS 2000 in Proj4 definitionraster::plot(rc_cob_rasterizacao, col = viridis::viridis(10))
plot(rc_cob$geom, add = TRUE)
Rasterização da cobertura da terra para Rio Claro/SP.
Vetorização
A operação inversa à rasterização é a vetorização, na qual iremos converter um raster em um vetor, sendo que esse vetor irá receber os valores dos pixels. O vetor em questão pode ser pontos (geralmente um gride de pontos), linhas (geralmente isolinhas ou linhas de contorno), ou polígonos (podendo esses polígonos ser ou não dissolvidos pelos valores dos pixels). Existem funções específicas para cada uma dessas conversões, sendo elas: raster::rasterToPoints(), raster::rasterToContour() e raster::rasterToPolygons(), respectivamente. Para a última função, ainda dispomos de uma alternativa mais veloz spex::polygonize().
Em nosso exemplo, vamos vetorizar o raster de elevação para Rio Claro/SP, criando um gride de pontos, sendo os pontos os centroides de cada pixels @ref(fig:fig-raster-vetorizacao-pontos)).
# vetorizacao de pontos
ra_rc_sirgas2000_utm23s_agre_media_pontos <- raster::rasterToPoints(ra_rc_sirgas2000_utm23s_agre_media, spatial = TRUE) %>%
sf::st_as_sf()raster::plot(ra_rc_sirgas2000_utm23s_agre_media, col = viridis::viridis(10, alpha = .8))
plot(ra_rc_sirgas2000_utm23s_agre_media_pontos, pch = 20, cex = .7, main = FALSE, add = TRUE)
plot(rc_2019_sirgas2000_utm23s$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Vetorização do raster de elevação criando um gride de pontos para Rio Claro/SP.
Nesse outro exemplo, vamos vetorizar o raster de elevação para Rio Claro/SP novamente, mas agora criando isolinhas @ref(fig:fig-raster-vetorizacao-linhas)).
# vetorizacao de linhas
ra_rc_sirgas2000_utm23s_agre_media_linhas <- raster::rasterToContour(x = ra_rc_sirgas2000_utm23s_agre_media) %>%
sf::st_as_sf()raster::plot(ra_rc_sirgas2000_utm23s_agre_media, col = viridis::viridis(10, alpha = .8))
contour(ra_rc_sirgas2000_utm23s_agre_media, labcex = 1, main = FALSE, add = TRUE)
plot(rc_2019_sirgas2000_utm23s$geom, col = NA, border = "red", lwd = 2, add = TRUE)
Vetorização do raster de elevação criando isolinhas para Rio Claro/SP.
Por fim, vamos vetorizar o raster de cobertura da terra criado anteriormente para Rio Claro/SP, criando polígonos não dissolvendo e dissolvidos @ref(fig:fig-raster-vetorizacao-poligonos)).
# vetorizacao de poligonos
rc_cob_rasterizacao_poligonos <- raster::rasterToPolygons(rc_cob_rasterizacao) %>%
sf::st_as_sf()
# vetorizacao de poligonos dissolvendo
rc_cob_rasterizacao_poligonos_dissolvidos <- raster::rasterToPolygons(rc_cob_rasterizacao, dissolve = TRUE) %>%
sf::st_as_sf()par(mfrow = c(1, 2))
raster::plot(rc_cob_rasterizacao, col = viridis::viridis(10))
plot(rc_cob_rasterizacao_poligonos$geometry, col = NA, border = "gray", lwd = 1, main = FALSE, add = TRUE)
raster::plot(rc_cob_rasterizacao, col = viridis::viridis(10))
plot(rc_cob_rasterizacao_poligonos_dissolvidos$geometry, col = NA, border = "gray", lwd = 1, main = FALSE, add = TRUE)
Vetorização do raster de cobertura da terra para Rio Claro/SP, não dissolvendo e dissolvendos os polígonos gerados.
par(mfrow = c(1, 1))Visualização de dados geográficos - Mapas
Um dos pontos finais de toda a análise envolvendo a manipulação de dados geográficos será a apresentação de um mapa com as informações de interesse espacializadas. Mas antes, é necessário ter conhecimento de alguns dos elementos principais para a composição de um mapa relativamente bem informativo. Além disso, o R nos permite criar tipos diferentes de mapas: estáticos, animados e interativos. Os mais comuns são os estáticos, mas podemos por vezes melhorar a apresentação dos dados geográficos criando mapas animados e/ou interativos, com o auxílio de páginas web. Por fim, veremos as melhores formas de exportar mapas para diferentes formatos.
Principais elementos de um mapa
Um mapa pode ser composto de vários elementos, tendo estes os intuito de auxiliar a visualização e entendimento de seu conteúdo. Apesar disso, nem todos os elementos necessitam estar presentes em todos os layouts de mapas, sendo que os mesmos devem atendem à necessidade das representações, podendo ser muitas vezes omitidos.
Os principais elementos de um mapa geralmente são compostos por: 1) mapa principal (ocupando quase toda a área da figura), 2) mapa secundário (geralmente muito menor que o mapa principal e com o intuito de mostrar a localização do mapa principal num contexto mais amplo, como país ou continente), 3) título (para resumir o intuito do mapa), 4) legenda (apresentando as informações detalhadas das classes ou escala de valores, geralmente identificando as cores e/ou texturas), barra de escala (representando quantas unidades do mapa representam do mundo real), 5) indicador de orientação (Norte) (indicando o norte geográfico, podendo ser representado por uma flecha, bússula ou compasso), 6) gride de coordenadas (coordenadas presentes nas laterais), 7) descrição do CRS (indicando qual o CRS), 8) informações de origem (informações sobre a fonte dos dados representados no mapa), 9) além de outros elementos auxiliares (como elementos textuais e figuras extras).
Podemos visualizar todos esses elementos resumidos na Figura @ref(fig:fig-mapa-elementos).
## Using year 2019## Scale bar set for latitude km and will be different at the top and bottom of the map.
Principais elementos de um mapa.
Principais pacotes para a composição de mapas
Há uma grande quantidade de pacotes para a composição de mapas no R. Aqui listamos os principais (Tabela @ref(tab:tab-mapa-pacotes)).
| Pacote | Descrição |
|---|---|
| ggplot2 | Cria visualizações de dados elegantes usando a gramática de gráficos |
| ggspatial | Estrutura de dados espaciais para ggplot2 |
| ggmap | Visualização espacial com ggplot2 |
| tmap | Mapas temáticos |
| leaflet | Cria mapas da web interativos com a biblioteca JavaScript ‘Leaflet’ |
| plotly | Cria gráficos interativos da Web por meio de ‘plotly.js’ |
| cartography | Cartografia temática |
| googleway | Cartografia temática |
| mapview | Acessa APIs do Google Maps para recuperar dados e mapas de plotagem |
| rasterVis | Visualização interativa de dados espaciais em R |
| cartogram | Métodos de visualização para dados raster |
| mapsf | Crie cartogramas com R |
| geogrid | Transforme polígonos geoespaciais em grades regulares ou hexagonais |
| geofacet | ‘ggplot2’ Utilitários de facetação para dados geográficos |
| globe | Plot 2D and 3D Views of the Earth, Including Major Coastline |
| linemap | Line Maps |
Mapas estáticos
Mapas estáticos são mapas simples e fixos para visualização de dados, sendo o tipo mais comum de saída visual. No início da composição de mapas no R, esse era o único tipo de mapa que a linguagem permitia produzir, principalmente utilizando o pacote sp (Pebesma e Bivand 2005). No entanto, com o advento de ferramentas de visualização dinâmicas no R, como componentes HTML, os mapas puderam ser compostos de forma dinâmica (animados e interativos).
Neste tópico abordaremos funções simples para composição de mapas estáticos, como o plot(), além de pacotes para composição de mapas mais elaborados, como os pacotes ggplot2 (Wickham 2016) e tmap (Tennekes 2018).
Função plot()
A função genérica plot() é a maneira mais rápida de compor mapas estáticos utilizando objetos espaciais vetoriais e raster, funcionando para ambos os pacotes que apresentamos anteriormente (sf e raster). Apesar da simplicidade, essa função geralmente tende a criar mapas com relativa velocidade, nos auxiliando principalmente em fases iniciais de desenvolvimento de um projeto. Essa função oferece dezenas de argumentos em base R, permitindo alguns ajustes limitados, com resultados bastante interessantes.
Como dito anteriormente, a função plot() vai funcionar diferentemente dependendo da classe do objeto espacial. Para objetos espaciais sf, a função vai plotar um mapa para cada coluna da tabela de atributos. Vamos usar de exemplo nosso mapa de biomas mostrado com os principais elementos de um mapa, podendo inclusive selecionar apenas a coluna de características geoespaciais (geom).
Primeiramente, vamos fazer o download dos dados de limites de biomas, retirando os sistemas costeiros, usando o pacote geobr ().
# biomas
biomas <- geobr::read_biomes(showProgress = FALSE) %>%
dplyr::filter(name_biome != "Sistema Costeiro")## Using year 2019Agora, quando utilizamos a função plot() para um objeto da classe sf, temos os três mapas, cada um indicando uma coluna da tabela de atritos (Figura @ref(fig:fig-vetor-map-plot)).
plot(biomas)
Mapa feito com a função plot() de um objeto sf para os Biomas do Brasil.
Selecionando as colunas desse objeto, podemos escolher a informação que queremos plotar, por exemplo, apenas a geometria geom. Além disso, podemos acrescentar os argumentos col para colorir e main para o título, além dos argumentos axes e graticule para adicionar as coordenadas e quadrículas, repectivamente. A legenda pode ser adicionada com a função legend() (Figura @ref(fig:fig-vetor-map-plot-selecao)).
plot(biomas$geom,
col = c("darkgreen", "orange", "orange4", "forestgreen", "yellow", "yellow3"),
main = "Biomas do Brasil", axes = TRUE, graticule = TRUE)
legend(x = -75, y = -20, pch = 15, cex = .7, pt.cex = 2.5, legend = biomas$name_biome,
col = c("darkgreen", "orange", "orange4", "forestgreen", "yellow", "yellow3"))
Para a classe dos objetos espaciais raster, a função plot() vai plotar um mapa para o tipo RasterLayer e quantos mapas houverem no objeto e couberem no espaço de plot para RasterBrick e RasterStack. Além disso, para essas classes do objeto raster, essa função provê também uma legenda e uma escala de cores automática (“terrain”). Vamos fazer o mapa da camada raster de elevação para o limite do município de Rio Claro/SP (Figura @ref(fig:fig-raster-layer-mapa)).
plot(ra_rc)
plot(rc_2019[1], col = NA, lwd = 2, add = TRUE)
Mapa feito com a função plot() de um objeto raster com uma camada.
Agora vamos plotar objetos da classe RasterStack, alterando a cor para “viridis”, usando a função viridis::viridis() do pacote homônimo. Vamos fazer o mapa de duas camadas raster bioclimáticas para o mundo (Figura @ref(fig:fig-raster-stack-mapa)).
plot(st[[1:2]], col = viridis::viridis(10))
Mapa feito com a função plot() de um objeto raster com várias camadas.
Para exportar esses mapas podemos utilizar as funções png() ou pdf(), indicando os argumentos para ter as configurações que desejamos, e finalizando com a função dev.off(). Vamos exportar, à título de exemplo, a última figura.
# diretorio
dir.create(here::here("dados", "mapas"))
# exportar mapa
png(filename = here::here("dados", "mapas", "elev_rc.png"), width = 20, height = 20, units = "cm", res = 300)
plot(ra_rc)
plot(rc_2019[1], col = NA, lwd = 2, add = TRUE)
dev.off()Pacotes ggplot2 e ggspatial
Como discutimos no capítulo xx sobre gráficos, o pacote ggplot2 utiliza a gramática de gráficos para composição de figuras no R (Wilkinson 1999, Wickhan 2016). Para cada classe de objeto geográfico há funções específicas para os dados: para objetos sf geom_sf() e para objetos raster geom_raster().
Além do pacote ggplot2, podemos utilizar o pacote ggspatial para acrescentar elementos geográficos como a barra de escala e o indicador de orientação (Norte), através das funções annotation_scale() e annotation_north_arrow(), respectivamente, além de outras funções específicas que não abordaremos nesta seção.
A estrutura de composição das funções do pacote ggplot2 vai funcionar parecido com a estruturação de gráficos já vista no capítulo XX, de modo que a cada função iremos utilizando o sinal de + para acrescentar outra camada. Iremos indicar os dados com a função ggplot() e a coluna da tabela de atributos que queremos representar com a função aes(). Em seguida, utilizamos a função geom_sf() para indicar que trata-se de um objeto sf.
Além dessas funções, podemos ainda fazer alterações nos mapas através das funções: scale_*() que vai alterar as características indicadas em aes(), coord_*() que vai alterar construção do mapa em relação às coordenadas, facet_*() que altera a disposição de vários mapas, e theme_*() e theme() que irão alterar características relacionadas ao tema, como cores de fundo, fontes e legenda. Podemos ainda utilizar as funções annotate() para adicionar textos e labs() para alterar o texto do título, legenda e eixos.
Vamos demonstrar esse funcionamento para compor o mapa de biomas, apresentado no início desta seção (Figura @ref(fig:fig-vetor-mapa-ggplot2)).
map_biomas_ggplot2 <- ggplot(data = bi) +
aes(fill = name_biome) +
geom_sf(color = "black") +
scale_fill_manual(values = c("darkgreen", "orange", "orange4",
"forestgreen", "yellow", "yellow3")) +
annotation_scale(location = "br") +
annotation_north_arrow(location = "br", which_north = "true",
pad_x = unit(0, "cm"), pad_y = unit(.5, "cm"),
style = north_arrow_fancy_orienteering) +
annotate(geom = "text", label = "CRS: SIRGAS2000/Geo", x = -38, y = -31, size = 2.5) +
annotate(geom = "text", label = "Fonte: IBGE (2019)", x = -39, y = -32.5, size = 2.5) +
labs(title = "Biomas do Brasil", fill = "Legenda", x = "Longitude", y = "Latitude") +
theme_bw() +
theme(title = element_text(size = 15, face = "bold"),
legend.title = element_text(size = 10, face = "bold"),
legend.position = c(.15, .25),
legend.background = element_rect(colour = "black"),
axis.title = element_text(size = 10, face = "plain"))
map_biomas_ggplot2## Scale on map varies by more than 10%, scale bar may be inaccurate
Mapa de Biomas do Brasil com o pacote ggplot2.
Para objetos raster, o uso do pacote ggplot2 para compor mapas requer um passo preliminar. Primeiramente, vamos criar um dataframe com os dados do raster, com as linhas sendo os pixels e as colunas sendo as coordenadas centrais da longitude e latitude, além dos valores de cada camada em cada coluna. Esse passo pode ser realizado com a função raster::rasterToPoints().
Uma vez que temos esses dados organizados, podemos utilizar as funções ggplot() para indicar o dataframe, e as colunas com a função aes(). Em seguida, utilizamos a função geom_raster() para indicar que trata-se de um objeto raster. Além dessas funções, podemos ainda utilizar as demais funções para alterar as características do mapa, como comentamos acima. Entretanto, devemos nos atentar para a função coord_*() e escolher aquela que vai fazer a construção do mapa em relação à coordenadas e resolução das células.
Como exemplo, vamos compor o mapa de elevação para Rio Claro/SP, adicionando também o limite do município (Figura @ref(fig:fig-raster-mapa-ggplot2)).
# dados
ra_rc_da <- raster::rasterToPoints(ra_rc) %>%
tibble::as_tibble()
head(ra_rc_da)## # A tibble: 6 x 3
## x y elevacao
## <dbl> <dbl> <dbl>
## 1 -47.8 -22.2 859
## 2 -47.8 -22.2 856
## 3 -47.8 -22.2 856
## 4 -47.8 -22.2 856
## 5 -47.8 -22.2 853
## 6 -47.8 -22.2 852# mapa
map_elev_rc_ggplot2 <- ggplot() +
geom_raster(data = ra_rc_da, aes(x = x, y = y, fill = elevacao)) +
geom_sf(data = rc_2019, color = "red", fill = NA, size = 1.3) +
scale_fill_viridis_c() +
coord_sf() +
annotation_scale(location = "br",
pad_x = unit(.5, "cm"), pad_y = unit(.7, "cm"),) +
annotation_north_arrow(location = "br", which_north = "true",
pad_x = unit(.4, "cm"), pad_y = unit(1.3, "cm"),
style = north_arrow_fancy_orienteering) +
annotate(geom = "text", label = "CRS: WGS84/Geo", x = -47.51, y = -22.53, size = 3) +
labs(title = "Elevação de Rio Claro/SP", fill = "Elevação (m)", x = "Longitude", y = "Latitude") +
theme_bw() +
theme(title = element_text(size = 15, face = "bold"),
legend.title = element_text(size = 10, face = "bold"),
legend.position = c(.2, .25),
legend.background = element_rect(colour = "black"),
axis.title = element_text(size = 10, face = "plain"),
axis.text.y = element_text(angle = 90, hjust = .4))
map_elev_rc_ggplot2
Mapa raster com o pacote ggplot2.
Para exportar mapas criados com o pacote ggplot2, podemos utilizar a função ggplot2::ggsave(), indicando os argumentos para ter as configurações que desejamos. Vamos exportar, à título de exemplo, a última figura.
# exportar mapa ggplot2
ggsave(filename = here::here("dados", "mapas", "elev_rc_ggplot2.png"), plot = map_elev_rc_ggplot2, width = 20, height = 20, units = "cm", dpi = 300)Pacote tmap
O pacote tmap é um pacote direcionado a criação de mapas temáticos, com uma sintaxe concisa que permite a criação de mapas com o mínimo de códigos, mas muito similar à sintaxe do pacote ggplot2 (Tennekes 2018). Ele também pode gerar mapas estáticos ou interativos usando o mesmo código, apenas mudando a forma de visualização com a função tmap_mode(), com o argumento mode igual a “plot” para estático e “view” para interativo. Por fim, o pacote tmap aceita diversas classes espaciais, incluindo objetos raster, de forma bastante simples. Mais sobre o pacote pode ser lido aqui. Novamente, atentar para a instalação extra em Linux e MacOS.
Todas as funções do pacote tmap iniciam-se com tm_*, facilitando seu uso. A cada função iremos utiliz o sinal de + para acrescentar outra camada, da mesma forma que o pacote ggplot2. A principal função, em que todos os objetos espaciais são dados de entrada, é tm_shape(). A partir dela, podemos seguir com funções específicas para vizualização de objetos sf, como tm_polygons(), tm_borders(), tm_fill(), tm_lines(), tm_dots() ou tm_bubbles(); ou com funções para objetos raster como tm_raster(). Ainda há funções como tm_text() para representação de textos das colunas da tabela de atributos, e tm_scale_bar(), tm_compass() e tm_graticules(), para adicionar barra de escala, indicador de orientação (Norte) e gride de coordenadas, repectivamente. Por fim, a função tm_credits() adiciona um texto descritivo e a função tm_layout() faz diversas mudanças nos detalhes e apresentação do mapa.
Uma funcionalidade muito interessante do pacote tmap é o uso da função tmaptools::palette_explorer() para escolher as paletas de cores disponíveis. Essa função requer que os pacotes shiny e shinyjs estejam instalados, e quando executada, retorna uma aba onde é possível editar e escolher algumas paletas de cores nativas do tmap.
Diversos parâmetros podem ser acrescentados às funções de composição do tmap, mas não as detalharemos aqui, pois todas são descritas nos vignettes do pacote: tmap: get started! e tmap: version changes.
Vamos seguir com a composição do mapa de biomas para o Brasil apresentado no início dessa seção (Figura @ref(fig-vetor-mapa-tmap)).
map_biomas_tmap <- tm_shape(bi, bbox = c(-74, -35, -27, 10)) +
tm_polygons(col = "name_biome",
pal = c("darkgreen", "orange", "orange4", "forestgreen", "yellow", "yellow3"),
border.col = "black",
title = "Legenda") +
tm_compass() +
tm_scale_bar(text.size = .6) +
tm_graticules(lines = FALSE) +
tm_credits("CRS: SIRGAS2000/Geo", position = c(.63, .13)) +
tm_credits("Fonte: IBGE (2019)", position = c(.63, .09)) +
tm_layout(title = "Biomas do Brasil",
title.position = c(.25, .95),
title.size = 1.8,
title.fontface = "bold",
legend.frame = TRUE,
legend.position = c("left", "bottom"),
legend.title.fontface = "bold")
map_biomas_tmap## Scale bar set for latitude km and will be different at the top and bottom of the map.
Mapa de Biomas do Brasil com o pacote tmap.
Além disso, o pacote tmap nos permite adicionar de forma simples um mapa secundário, provendo uma localização regional de interesse (Figura @ref(fig:fig-vetor-mapa-sec-tmap)).
# dados
sa <- rnaturalearth::ne_countries(continent = "South America")
br <- rnaturalearth::ne_countries(country = "Brazil")
bi <- geobr::read_biomes(showProgress = FALSE) %>%
dplyr::filter(name_biome != "Sistema Costeiro")## Using year 2019# mapa secundario
map_biomas_sa <- tm_shape(sa) +
tm_polygons() +
tm_shape(br) +
tm_polygons(col = "gray50")
# juntando os mapas
map_biomas_tmap## Scale bar set for latitude km and will be different at the top and bottom of the map.print(map_biomas_sa, vp = viewport(.815, .875, wi = .2, he = .2))Mapa vetorial primário e secundário com o pacote tmap.
Como exemplo de mapa raster, vamos compor novamente o mapa de elevação para Rio Claro/SP, adicionando também o limite do município (Figura @ref(fig:fig-raster-mapa-tmap)).
map_elev_rc_tmap <- tm_shape(ra_rc) +
tm_raster(pal = "viridis", title = "Elevação (m)") +
tm_shape(rc_2019) +
tm_borders(col = "red", lwd = 2) +
tm_compass(position = c(.9, .08)) +
tm_scale_bar(text.size = .6, position = c(.67, 0)) +
tm_graticules(lines = FALSE) +
tm_credits("CRS: WGS84/Geo", position = c(.67, .06)) +
tm_layout(title = "Elevação Rio Claro/SP",
title.size = 1.3,
title.fontface = "bold",
legend.title.size = .7,
legend.text.size = .6,
legend.frame = TRUE,
legend.position = c(.01, .01),
legend.title.fontface = "bold")
map_elev_rc_tmap
Mapa raster de elevação com o pacote tmap.
Para exportar mapas criados com o pacote tmap podemos utilizar a função tmap::tmap_save(), indicando os argumentos para ter as configurações que desejamos. Vamos exportar, à título de exemplo, a última figura.
# exportar mapa tmap
tmap::tmap_save(tm = map_elev_rc_tmap, filename = here::here("dados", "mapas", "elev_rc_tmap.png"), width = 20, height = 20, units = "cm", dpi = 300)Mapas animados
Podemos montar mapas facetados para mostrar como padrões espaciais variam ao longo do tempo, como por exemplo, os limites do Brasil ao longo dos anos (Figura @ref(fig:fig-vetor-brasil-tmap-facted)). Entretanto essa a abordagem possui algumas desvantagens, de modo que as facetas podem ficar muito pequenas quando há muitas delas.
# dados
br_anos <- NULL
for(i in c(1872, 1900, 1911, 1920, 1933, 1940, 1950, 1960, 1970, 1980, 1991, 2001, 2010, 2019)){
br_anos <- geobr::read_state(code_state = "all", year = i, showProgress = FALSE) %>%
dplyr::mutate(year = i) %>%
dplyr::bind_rows(br_anos, .)
}## Using year 1872## Loading data for the whole country## Using year 1900## Loading data for the whole country## Using year 1911## Loading data for the whole country## Using year 1920## Loading data for the whole country## Using year 1933## Loading data for the whole country## Using year 1940## Loading data for the whole country## Using year 1950## Loading data for the whole country## Using year 1960## Loading data for the whole country## Using year 1970## Loading data for the whole country## Using year 1980## Loading data for the whole country## Using year 1991## Loading data for the whole country## Using year 2001## Loading data for the whole country## Using year 2010## Loading data for the whole country## Using year 2019## Loading data for the whole country# numero de estados ao longo do tempo
br_anos$year %>% table## .
## 1872 1900 1911 1920 1933 1940 1950 1960 1970 1980 1991 2001 2010 2019
## 21 21 22 22 22 24 28 29 28 28 29 27 27 27# mapa facetado
map_brasil_tmap <- tm_shape(br_anos) +
tm_polygons() +
tm_facets(by = "year", nrow = 4)
map_brasil_tmap
Mapa vetor facetado dos estados brasileiros ao longo do tempo com o pacote tmap.
Uma solução é a composição de mapas animados. Apesar de dependerem da publicação digital, os mapas animados podem aprimorar relatórios físicos à medida que o vínculo a uma página da web contendo a versão animada torna-se simples. Existem várias maneiras de gerar animações em R, e uma forma é com o pacote gganimate e ggplot2. Entretanto, aqui veremos a criação de mapas animados com tmap.
Podemos criar mapas animados alterando dois argumentos da função tm_facets():
- trocando o by = “year” por along = “year”
- indicando o free.coords = FALSE
Por fim, podemos exportar o mapa animado no formato de .gif utilizando a função tmap::tmap_animation(), indicando a taxa de atualização com o argumento delay (Figura @ref(fig:fig-vetor-brasil-tmap-animated)). Alguns pacotes extras são requeridos dependendo do sistema operacional utilizado.
# mapa animado
map_brasil_tmap_ani <- tm_shape(br_anos) +
tm_polygons() +
tm_facets(along = "year", free.coords = FALSE)
# exportar
tmap::tmap_animation(tm = map_brasil_tmap_ani, filename = here::here("dados", "mapas", "elev_rc_tmap_ani.gif"), delay = 30)
Mapa vetorial animado mostrando os estados brasileiros ao longo do tempo com o pacote tmap.
Mapas interativos
Mapas interativos podem assumir muitas formas, sendo que a mais comum e útil é a capacidade de deslocar e ampliar qualquer parte de um conjunto de dados geográficos sobreposto em um “mapa da web”. Diversos pacotes nos permitem criar esse tipo de mapa, sendo os mais comuns o tmap, mapview e leaflet. É importante destacar ainda que esses mapas irão ser compostos numa janela especial de “Viewer”.
Pacote tmap
Um recurso exclusivo do tmap é sua capacidade de criar mapas estáticos e interativos usando o mesmo código. Os mapas podem ser visualizados interativamente em qualquer ponto mudando para o modo de visualização, usando a função tmap::tmap_mode(mode = "view") (Figura @ref(fig:fig-raster-mapa-tmap-int)).
# mudar o modo de exibicao do tmap
tmap::tmap_mode(mode = "view")
map_elev_rc_tmap_int <- map_elev_rc_tmap
map_elev_rc_tmap_intMapa vetorial interativo com o pacote tmap.
Para exportar mapas interativos criados com o pacote tmap, podemos utilizar novamente a função tmap::tmap_save(), indicando a extensão como .html.
# exportar mapa tmap interativo
tmap::tmap_save(tm = map_elev_rc_tmap_int,
filename = here::here("dados", "mapas", "elev_rc_tmap_int.html"))Pacote mapview
O pacote mapview cria rapidamente mapas interativos simples com a função mapvew::mapview() (Figura @ref(fig:fig-raster-mapa-mapview-int)). Entretanto, outras características podem ser mudadas para criar mapas bem mais elaborados, como pode ser visto através do site do pacote.
map_elev_rc_mapview_int <- mapview::mapview(ra_rc, col.regions = viridis::viridis(100))
map_elev_rc_mapview_intMapa vetorial interativo com o pacote mapview.
Para exportar mapas interativos criados com o pacote mapview, podemos utilizar a função mapivew::mapshot(), indicando a extensão como .html.
# exportar mapa tmap interativo
mapview::mapshot(x = map_elev_rc_mapview_int,
url = here::here("dados", "mapas", "elev_rc_mapview_int.html"))Pacote leaflet
O leaflet é o pacote de mapeamento interativo mais utilizado e completo em R. Esse pacote fornece uma interface utilizando a biblioteca JavaScript e muitos argumentos podem ser compreendidos lendo a documentação da biblioteca original.
Mapas interativos usando esse pacote são criados utilizando a função leaflet::leaflet(). O resultado dessa função é um objeto da classe leaflet, que pode ser alterado por outras funções deste pacote, permitindo que várias camadas e configurações de controle sejam adicionadas interativamente (Figura @ref(fig:fig-raster-mapa-leaflet-int)).
Mais sobre o pacote leaflet pode ser consultado em seu site e CheatSheet.
# paleta de cores
pal <- colorNumeric(viridis::viridis(10), raster::values(ra_rc))
# mapa
map_elev_rc_leaflet_int <- leaflet() %>%
addProviderTiles("CartoDB.Positron") %>%
addRasterImage(ra_rc, colors = pal, opacity = .8) %>%
addLegend(pal = pal, values = raster::values(ra_rc), title = "Elevação (m)") %>%
addPolygons(data = rc_2019, col = "red", fill = NA)
map_elev_rc_leaflet_intMapa vetorial interativo com o pacote leaflet.
Para exportar mapas interativos criados com o pacote leaflet, podemos utilizar novamente a função mapivew::mapshot(), indicando a extensão como .html.
# exportar mapa tmap interativo
mapview::mapshot(x = map_elev_rc_leaflet_int,
url = here::here("dados", "mapas", "elev_rc_leaflet_int.html"))Exemplos de aplicações de análises geográficas para dados ecológicos
Agora que vimos os conceitos e aplicações básicas de manejo e visualização de dados geográficos, podemos avançar para realizar três exemplos de aplicações para dados ecológicos. Para isso, usaremos novamente os dados de comunidades de anfíbios da Mata Atlântica (Atlantic Amphibians, Vancine et al. 2018). Primeiramente, veremos como resumir informações de biodiversidade (número de ocorrências e riqueza) para hexágonos. Num segundo momento, veremos como associar dados ambientais para coordenadas de espécies ou comunidades. Por fim, iremos realizar predições espaciais contínuas de adequabilidade de habitat e número de espécies.
Resumir informações de biodiversidade para unidades espaciais
Resumir informações para unidades espaciais é um passo muito frequente em análises de Macroecologia, Biogeogradia ou Ecologia da Paisagem. Nesta seção, iremos contabilizar o número de ocorrências e riqueza de anfíbios para hexágonos na Mata Atlântica.
Primeiramente, vamos importar e preparar os dados de biodiversidade que usaremos nesses exemplos. Vamos começar importanto os locais de amostragens de anfíbios na Mata Atlântica e selecionando apenas as colunas de interesse.
# importar locais
aa_locais <- readr::read_csv(here::here("dados", "tabelas", "ATLANTIC_AMPHIBIANS_sites.csv")) %>%
dplyr::select(id, longitude, latitude, species_number)##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_character(),
## reference_number = col_double(),
## species_number = col_double(),
## month_start = col_double(),
## year_start = col_double(),
## month_finish = col_double(),
## year_finish = col_double(),
## effort_months = col_double(),
## latitude = col_double(),
## longitude = col_double(),
## altitude = col_double(),
## temperature = col_double(),
## precipitation = col_double()
## )
## ℹ Use `spec()` for the full column specifications.Agora vamos importar as espécies das comunidades, selecionando apenas as espécies com nomes válidos, e transformando a coluna de indivíduos para 1, para compor posteriormente uma matriz de comunidade de espécies.
# importar especies
aa_especies <- readr::read_csv(here::here("dados", "tabelas", "ATLANTIC_AMPHIBIANS_species.csv")) %>%
tidyr::drop_na(valid_name) %>%
dplyr::select(id, valid_name, individuals) %>%
dplyr::distinct(id, valid_name, .keep_all = TRUE) %>%
dplyr::mutate(individuals = tidyr::replace_na(individuals, 1),
individuals = ifelse(individuals > 0, 1, 1))##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## id = col_character(),
## order = col_character(),
## superfamily = col_character(),
## family = col_character(),
## subfamily = col_character(),
## genus = col_character(),
## species = col_character(),
## valid_name = col_character(),
## individuals = col_double(),
## endemic = col_double()
## )Podemos agora juntar a tabela de locais, que possui as coordenadas à tabela de espécies. Em seguida convertemos essa única tabela na classe vetor sf.
# join das coordenadas e classe sf
aa_especies_locais_ve <- aa_especies %>%
dplyr::left_join(aa_locais) %>%
dplyr::relocate(longitude, latitude, .after = 1) %>%
dplyr::mutate(lon = longitude, lat = latitude) %>%
sf::st_as_sf(coords = c("lon", "lat"), crs = 4326)## Joining, by = "id"Agora vamos baixar o limite do Bioma da Mata Atlântica para o Brasil, converter o GCS para WGS84/Geo e ajustar sua extensão para remover as ilhas no Oceano Atlântico.
# mata atlantica
mata_atlantica <- geobr::read_biomes(showProgress = FALSE) %>%
dplyr::filter(name_biome == "Mata Atlântica") %>%
sf::st_transform(crs = 4326) %>%
sf::st_crop(xmin = -55, ymin = -30, xmax = -34, ymax = -5)## Using year 2019## although coordinates are longitude/latitude, st_intersection assumes that they are planar## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriesPodemos vericar se as coordenadas e o limite do bioma estão todos corretos compondo um mapa preliminar, usando o pacote tmap (Figura @ref(fig:fig-vetor-aa-ma)).
# mapa
tm_shape(mata_atlantica, bbox = aa_especies_locais_ve) +
tm_polygons() +
tm_shape(aa_especies_locais_ve) +
tm_dots(size = .1, col = "forestgreen")
Mapa dos locais do Atlantic Amphibians e do limite da Mata Atlântica.
Como o limite utilizado para reunir informações das comunidade de anfíbios foi o mais abrangente possível (Muylaert et al. 2018, Vancine et al. 2018), iremos selecionar apenas os locais que caem dentro do limite da Mata Atlântica que estamos utilizando aqui.
# selecionar os locais dentro do limite
aa_especies_locais_ve_ma <- aa_especies_locais_ve[mata_atlantica, ]## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planarPodemos refazer o mapa mostrando as coordenadas retiradas em vermelho e as que ficaram em verde (Figura @ref(fig:fig-vetor-aa-ma-sel)).
# mapa
tm_shape(mata_atlantica, bbox = aa_especies_locais_ve) +
tm_polygons() +
tm_shape(aa_especies_locais_ve) +
tm_bubbles(size = .1, col = "red") +
tm_shape(aa_especies_locais_ve_ma) +
tm_bubbles(size = .1, col = "forestgreen")
Mapa dos locais do Atlantic Amphibians que caem dentro do limite da Mata Atlântica.
O próximo passo é criar um gride de hexágonos para o Bioma da Mata Atlântica. Usaremos a função sf::st_make_grid() que pode criar quadrículas ou hexágonos. Esses hexágonos terão a área equivalente à quadríulas de 1º de tamanho (aproximadamente 110 km²). Usaremos a função sf::st_area() para calcular as áreas dos hexágonos e a função tibble::rowid_to_column() para criar uma identificação para cada feição.
# criar os hexagonos
mata_atlantica_hex <- sf::st_make_grid(x = mata_atlantica,
cellsize = 1,
square = FALSE) %>%
sf::st_as_sf() %>%
dplyr::mutate(areakm2 = sf::st_area(.)/1e6) %>%
tibble::rowid_to_column("id_hex")
# selecionar os hexagonos para dentro do limite da mata atlantica
mata_atlantica_hex <- mata_atlantica_hex[mata_atlantica,]## although coordinates are longitude/latitude, st_intersects assumes that they are planarPodemos conferir os hexágonos criados fazendo um mapa preliminar (Figura @ref(fig:fig-vetor-aa-ma-sel-hex)).
# mapa
tm_shape(mata_atlantica, bbox = mata_atlantica_hex) +
tm_polygons() +
tm_shape(mata_atlantica_hex) +
tm_borders()
Mapa dos hexágonos para o limite da Mata Atlântica.
Podemos agora associar as espécies aos hexágonos fazendo um “join” espacial, utilizando a função sf::st_join().
mata_atlantica_hex_especies <- sf::st_join(x = mata_atlantica_hex,
y = aa_especies_locais_ve_ma,
left = TRUE)## although coordinates are longitude/latitude, st_intersects assumes that they are planarPor fim, podemos agregar os dados para ter o número de ocorrências e de espécies por hexágono.
mata_atlantica_hex_especies_oco_riq <- mata_atlantica_hex_especies %>%
dplyr::group_by(id_hex) %>%
dplyr::summarise(ocorrencias = length(valid_name[!is.na(valid_name)]),
riqueza = n_distinct(valid_name, na.rm = TRUE))## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planarFinalmente podemos compor os mapas finais, mostrando os hexágonos com cores e valores do número de ocorrências e de espécies @ref(fig:fig-vetor-aa-ma-sel-hex-oco-riq)).
# mapa de ocorrencias
mapa_oco <- tm_shape(mata_atlantica_hex_especies_oco_riq) +
tm_polygons(title = "Ocorrência de anfíbios", col = "ocorrencias",
pal = "-RdYlBu", style = "pretty") +
tm_text("ocorrencias", size = .4) +
tm_graticules(lines = FALSE) +
tm_compass() +
tm_scale_bar() +
tm_layout(legend.title.size = 2,
legend.title.fontface = "bold",
legend.position = c("left", "top"))
# mapa de riqueza
mapa_riq <- tm_shape(mata_atlantica_hex_especies_oco_riq) +
tm_polygons(title = "Riqueza de anfíbios", col = "riqueza",
pal = "-Spectral", style = "pretty") +
tm_text("riqueza", size = .4) +
tm_graticules(lines = FALSE) +
tm_compass() +
tm_scale_bar() +
tm_layout(legend.title.size = 2,
legend.title.fontface = "bold",
legend.position = c("left", "top"))
# uniao dos mapas
tmap_arrange(mapa_oco, mapa_riq)## Scale bar set for latitude km and will be different at the top and bottom of the map.
## Scale bar set for latitude km and will be different at the top and bottom of the map.
Mapa com o número de ocorrências e riqueza de anfíbios para hexágonos no limite da Mata Atlântica.
Atribuição de dados climáticos para pontos
Atribuir informações ambientais à ocorrências é um passo fundamental para diversas análises. Nesta seção, iremos atribuir os valores das variáveis bioclimáticas aos locais de amostragem de anfíbios na Mata Atlântica.
Já realizamos o download das variáveis bioclimáticas na seção xx. Vamos importar novamente esses dados, primeiramente listando as camadas e depois importando com a função raster:stack().
# listar arquivos
fi <- dir(path = here::here("dados", "raster"), pattern = "wc") %>%
grep(".tif", ., value = TRUE)
# importar
var <- raster::stack(here::here("dados", "raster", fi))
# renomear
names(var) <- c("bio01", paste0("bio", 10:19), paste0("bio0", 2:9))Da seção antetior, já temos o objeto com a tabela de coordenadas dos locais de amostragem das comunidades de anfíbios. Vamos agora criar um objeto vetorial das coordenadas e em seguida selecionar os locais dentro do limite do bioma da Mata Atlântica.
# importar locais
aa_locais_ve <- aa_locais %>%
dplyr::mutate(lon = longitude, lat = latitude) %>%
sf::st_as_sf(coords = c("lon", "lat"), crs = 4326)
# selecionar locais para o limite
aa_locais_ve## Simple feature collection with 1163 features and 4 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -56.74194 ymin: -33.51083 xmax: -34.79667 ymax: -3.51525
## geographic CRS: WGS 84
## # A tibble: 1,163 x 5
## id longitude latitude species_number geometry
## * <chr> <dbl> <dbl> <dbl> <POINT [°]>
## 1 amp1001 -43.4 -8.68 19 (-43.42194 -8.68)
## 2 amp1002 -38.9 -3.55 16 (-38.85783 -3.545527)
## 3 amp1003 -38.9 -3.57 14 (-38.88869 -3.574194)
## 4 amp1004 -38.9 -3.52 13 (-38.9188 -3.51525)
## 5 amp1005 -38.9 -4.28 30 (-38.91083 -4.280556)
## 6 amp1006 -36.4 -9.23 42 (-36.42806 -9.229167)
## 7 amp1007 -40.9 -3.85 23 (-40.89444 -3.846111)
## 8 amp1008 -40.9 -3.83 19 (-40.91944 -3.825833)
## 9 amp1009 -40.9 -3.84 13 (-40.91028 -3.8375)
## 10 amp1010 -35.2 -6.14 1 (-35.22944 -6.136944)
## # … with 1,153 more rowsUsaremos agora a função raster::extract() para extrair e associar os valores das variáveis bioclimáticas para os locais de amostragem.
# extract
aa_locais_ve_var <- aa_locais_ve %>%
dplyr::mutate(raster::extract(var, ., df = TRUE)) %>%
dplyr::select(-ID) %>%
dplyr::relocate(bio02:bio09, .after = bio01)Podemos ver esses dados na Tabela @ref(tab:tab-aa-var).
| id | longitude | latitude | species_number | bio01 | bio02 | bio03 | bio04 | bio05 |
|---|---|---|---|---|---|---|---|---|
| amp1001 | -43.42194 | -8.680000 | 19 | 24.88622 | 12.842188 | 76.28720 | 84.54608 | 33.17875 |
| amp1002 | -38.85783 | -3.545527 | 16 | 26.43918 | 7.802419 | 78.70332 | 59.72862 | 31.37876 |
| amp1003 | -38.88869 | -3.574194 | 14 | 26.43918 | 7.802419 | 78.70332 | 59.72862 | 31.37876 |
| amp1004 | -38.91880 | -3.515250 | 13 | 26.43918 | 7.802419 | 78.70332 | 59.72862 | 31.37876 |
| amp1005 | -38.91083 | -4.280556 | 30 | 22.52411 | 8.434667 | 73.76507 | 64.62517 | 28.08925 |
| amp1006 | -36.42806 | -9.229167 | 42 | 21.61951 | 8.110271 | 67.74083 | 147.67953 | 27.98150 |
Podemos ainda fazer alguns mapas para espacializar essas variáveis (Figura @ref(fig:fig-vetor-aa-var)).
# mapa
aa_locais_ve_var %>%
dplyr::select(bio01:bio06) %>%
tidyr::gather(var, val, -geometry) %>%
tm_shape() +
tm_bubbles(size = .1, col = "val", pal = "-Spectral") +
tm_facets("var", free.scales = TRUE) +
tm_layout(legend.outside = FALSE)
Mapa mostrando os valores das variáveis bioclimáticas (BIO01:BIO06) para os locais amostrados de comunidades de anfíbios para hexágonos no limite da Mata Atlântica.
Predições espaciais de objetos raster
O pacote raster além de permitir realizar manejo e visualização de dados raster no R, também permite a extrapolação do ajuste de análises, como GLMs, GAMs dentre outras. Aqui, faremos uma pequena demostração utilizando a função raster::predict(), predizendo o resultado de dois ajustes de GLMs para a presença/ausência de uma espécie de anuro e a extrapolação do número de espécies de anfíbios para o Bioma da Mata Atlântica.
Para ajutar um GLM para dados de presença/ausência, podemos usar a tabela já criada anteriormente, com as espécies e as coordenadas, e fazer um join com a última tabela que criamos com os dados bioclimáticos.
aa_locais_ve_var_especies <- aa_especies %>%
dplyr::left_join(., sf::st_drop_geometry(aa_locais_ve_var), by = "id")Agora, vamos selecionar ocorrências da espécie Haddadus binotatus, atribuindo 1 quando ela ocorre e 0 quando ela não ocorre. Essa espécie é relativamente comum na serrapilheira de florestas da Mata Atlântica, e recebe esse nome em homenagem a um grande pesquisador de anfíbios da Mata Atlântica, o Prof. Célio Fernando Baptista Haddad.
aa_locais_ve_var_especies_hb <- aa_locais_ve_var_especies %>%
dplyr::mutate(pa = ifelse(valid_name == "Haddadus binotatus", 1, 0), .after = individuals) %>%
dplyr::distinct(id, .keep_all = TRUE)Vamos utilizar apenas as variáveis não correlacionadas para o índice de correlação de Pearson para r < 0,7.
# correlacao
corr <- aa_locais_ve_var_especies_hb %>%
dplyr::select(bio01:bio19) %>%
cor() %>%
caret::findCorrelation(.7, names = TRUE)
# selecao das variaveis nao correlacionadas
aa_locais_ve_var_especies_hb_cor <- aa_locais_ve_var_especies_hb %>%
dplyr::select(pa, bio01:bio19) %>%
dplyr::select(-c(corr))## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(corr)` instead of `corr` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.Agora sim, podemos ajustar um modelo simples da presença e ausência dessa espécie, utilizando as variáveis não correlacionadas, através de um GLM para a família binomial (para mais detalhes volte para o capítulo xx).
# ajustar glm
modelo_pa <- glm(formula = pa ~ ., data = aa_locais_ve_var_especies_hb_cor, family = binomial("logit"))## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredAntes de fazermos a predição da distribuição potencial da espécie é fundamental que o objeto raster esteja ajustado para o limite da Mata Atlântica. Para isso vamos utilizar as funções raster::crop() e raster::mask() para fazer esse ajuste @ref(fig:fig-raster-bio-ajuste)).
# ajuste da extensao e limite
var_mata_atlantica <- var %>%
raster::crop(mata_atlantica) %>%
raster::mask(mata_atlantica)# map
tm_shape(var_mata_atlantica[[c(1, 4)]]) +
tm_raster(pal = "viridis", title = c("bio01", "bio12")) +
tm_facets(free.scales.raster = TRUE)
Mapa de dois rasters (BIO01 e BIO12) ajustados ao limite da Mata Atlântica.
Agora podemos fazer a predição desse modelo para todo o bioma da Mata Atlântica. Essa função vai utilizar os coeficinetes do modelo ajustado para gerar um raster de predição para todos os pixels da Mata Atlântica.
# predicoes
modelo_pa_pred <- raster::predict(model = modelo_pa, object = var_mata_atlantica)
modelo_pa_pred## class : RasterLayer
## dimensions : 149, 121, 18029 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -55, -34.83333, -30, -5.166667 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : -27.40322, -1.396923 (min, max)Por fim, no último passo podemos usar as ocorrências para tornar esse modelo binário, ou seja, apenas com valores 0 ou 1. Para isso vamos extrair os valores da predição para as ocorrências da espécies e escolher o menor como sendo o corte. A partir desse valor consideraremos o pixels acima como 1 e abaixo como 0 (Pearson et al. 2007).
# menor valor da predicao
menor_valor <- modelo_pa_pred %>%
raster::extract(aa_locais_ve_var_especies_hb[, c("longitude", "latitude")]) %>%
min(na.rm = TRUE)
# selecao dos pixels de presenca/ausencia potencial
modelo_pa_pred_corte <- modelo_pa_pred >= menor_valorPor fim, vamos produzir dois mapas mostrando os valores das precições e o mapa binário. @ref(fig:fig-raster-pred-modelo-bin)).
# ocorrencias de hb
hb <- aa_locais_ve_var_especies_hb %>%
dplyr::filter(pa == 1) %>%
sf::st_as_sf(coords = c("longitude", "latitude"), crs = 4326)
# mapa predicao continua
mapa_pred_cont <- tm_shape(modelo_pa_pred) +
tm_raster(title = "Predição contínua", pal = "-Spectral") +
tm_shape(hb) +
tm_bubbles(size = .2) +
tm_graticules(lines = FALSE) +
tm_compass() +
tm_scale_bar() +
tm_layout(legend.title.size = 2,
legend.title.fontface = "bold",
legend.position = c("left", "top"))
# mapa predicao binaria
mapa_pred_bin <- tm_shape(modelo_pa_pred_corte) +
tm_raster(title = "Predição binária", pal = "-Spectral",
labels = c("Potencialmente ausente", "Potencialmente presente")) +
tm_shape(hb) +
tm_bubbles(size = .2) +
tm_graticules(lines = FALSE) +
tm_compass() +
tm_scale_bar() +
tm_layout(legend.title.size = 2,
legend.title.fontface = "bold",
legend.position = c("left", "top"))
# uniao dos mapas
tmap_arrange(mapa_pred_cont, mapa_pred_bin)## Scale bar set for latitude km and will be different at the top and bottom of the map.
## Scale bar set for latitude km and will be different at the top and bottom of the map.## Legend labels were too wide. The labels have been resized to 0.49, 0.48. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
Mapa da predição contínua e binária do modelo ajustado para a presença/ausência da espécie Haddadus binotatus na Mata Atlântica.
Em nossa segunda análise, vamos predizer os dados de riqueza para todo o bioma da Mata Atlântica. Para isso, temos de retirar novamente as variáveis correlacionadas.
# correlacao
corr <- aa_locais_ve_var %>%
sf::st_drop_geometry() %>%
dplyr::select(bio01:bio19) %>%
cor() %>%
caret::findCorrelation(.7, names = TRUE)
# selecao das variaveis nao correlacionadas
aa_locais_var_cor <- aa_locais_ve_var %>%
sf::st_drop_geometry() %>%
dplyr::select(species_number, bio01:bio19) %>%
dplyr::select(-c(corr))Agora sim, podemos criar os GLMs com famílias de distribuição apropriadas para dados de contabem como Poisson e Binomial Negativa.
# modelo poisson
modelo_riq_pois <- glm(formula = species_number ~ ., data = aa_locais_var_cor, family = poisson)
# modelo binomial negativo
modelo_riq_nb <- MASS::glm.nb(formula = species_number ~ ., data = aa_locais_var_cor)Com os modelos ajustados, podemos fazer as predições utilizando os objetos raster com as variáveis ambientais.
# predicao do modelo poisson
modelo_riq_pois <- predict(model = modelo_riq_pois, object = var_mata_atlantica)
modelo_riq_pois_contagem <- exp(modelo_riq_pois)
# predicao do modelo binomial negativo
modelo_riq_nb <- predict(model = modelo_riq_nb, object = var_mata_atlantica)
modelo_riq_nb_contagem <- exp(modelo_riq_nb)Por fim, podemos compor os dois mapas de predições (Figura @ref(fig:fig-raster-pred-modelo-riq)).
# mapa predicao poisson
mapa_pred_riq_pois <- tm_shape(modelo_riq_pois) +
tm_raster(title = "Número de espécies (Poisson)", pal = "-Spectral") +
tm_graticules(lines = FALSE) +
tm_compass() +
tm_scale_bar() +
tm_layout(legend.title.size = 2,
legend.title.fontface = "bold",
legend.position = c("left", "top"))
# mapa predicao binomial negativo
mapa_pred_riq_nb <- tm_shape(modelo_riq_nb) +
tm_raster(title = "Número de espécies (Binomial Negativa)", pal = "-Spectral") +
tm_graticules(lines = FALSE) +
tm_compass() +
tm_scale_bar() +
tm_layout(legend.title.size = 2,
legend.title.fontface = "bold",
legend.position = c("left", "top"))
# uniao dos mapas
tmap_arrange(mapa_pred_riq_pois, mapa_pred_riq_nb)## Scale bar set for latitude km and will be different at the top and bottom of the map.
## Scale bar set for latitude km and will be different at the top and bottom of the map.
Mapa da predição de riqueza utilizando o modelo Poisson e Binomial Negativa para a Mata Atlântica.
Para se aprofundar
Listamos aqui as principais referências sobre manipulação, visualização de dados geográficos e análises espaciais no R.
Lovelace, Nowosad & Muenchow (2019) Geocomputation with R. Chapman & Hall/CRC Biostatistics Series
Mas et al. (2019) Análise espacial com R.
Pebesma & Bivand (2020) Spatial Data Science.
Moraga, Paula (2019) Geospatial Health Data: Modeling and Visualization with R-INLA and Shiny. Chapman & Hall/CRC Biostatistics Series
Brunsdon & Comber (2019) An Introduction to Spatial Analysis and Mapping in R. 2nd edition.
Mieno (2020) R as GIS for Economists.
Gimond (2021) Intro to GIS and Spatial Analysis.
Engel (2019) Using Spatial Data with R.
Dorman (2021) Using R for Spatial Data Analysis
Rowe, Arribas-Bel (2021) Spatial Modelling for Data Scientists
Wegmann, Leutner & Dech (2016) Remote Sensing and GIS for Ecologists: Using Open Source Software
Wegmann, Schwalb-Willmann & Dech (2020) An Introduction to Spatial Data Analysis Remote Sensing and GIS with Open Source Software.
Fletcher & Fortin (2018) Spatial ecology and conservation modeling: Applications with r. Springer International Publishing.
Lapaine, Miljenko, Usery, E. Lynn (2017) Choosing a Map Projection. Spring.